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ABSTRACT 


Current  practices  for  modeling  the  ocean  floor  in 
underwater  explosion  simulations  call  for  application  of  an 
inviscid  fluid  with  soil  properties.  A  method  for  modeling 
the  ocean  floor  as  a  Lagrangian  solid,  vice  an  Eulerian 
fluid,  was  developed  in  order  to  determine  its  effects  on 
underwater  explosions  in  shallow  water  using  the  DYSMAS 
solver.  The  Lagrangian  solid  bottom  model  utilized 
transmitting  boundary  segments,  exterior  nodal  forces 
acting  as  constraints,  and  the  application  of  prestress  to 
minimize  any  distortions  into  the  fluid  domain.  Elastic 
materials  were  used,  though  multiple  constitutive  soil 
models  can  be  applied  to  improve  the  accuracy.  This  method 
is  unable  to  account  for  soil  cratering  effects,  however  it 
provides  the  distinct  advantage  of  modeling  contoured  ocean 
floors  such  as  dredged  channels  and  sloped  bottoms  absent 
in  Eulerian  formulations.  The  dynamic  loading  effects  of 
the  investigated  bottom  contours  were  found  to  be 
negligible  in  the  analyzed  cases  as  a  result  of  the  bulk 
cavitation  zone  which  dominates  the  chosen  fluid  field  and 
serves  as  a  buffer  to  the  target.  In  addition  to  its 
utility  in  bottom  modeling,  implementation  of  the  non¬ 
reflecting  boundary  along  with  realistic  material  models 
can  be  used  to  drastically  reduce  the  size  of  current  fluid 
domains . 
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I. 


INTRODUCTION 


A .  BACKGROUND 

Research  into  the  phenomena  of  underwater  explosions 
(UNDEX)  and  their  effects  on  ship  structure  and  equipment 
was  conducted  in  earnest  following  World  War  II.  The 
motivation  was  that  numerous  ships  experienced  extensive 
damage  not  only  from  direct  hits,  but  also  from  near  misses 
during  the  war.  The  foundation  of  the  field  of  study  is 
based  in  the  works  of  Arons  [1],  Cole  [2],  Snay  [3],  and 
Taylor  [4].  These  researchers  provided  extensive 
theoretical  and  empirical  relationships  describing  the 
various  phenomena  of  underwater  explosions.  Their  research 
led  to  the  development  of  shock  design  and  testing 
parameters  by  the  U.S.  Navy.  The  most  recent  set  of  these 
requirements  are  delineated  in  MIL-S-901D  [5],  NAVSEA  0908- 
LP-000-3010A  [6],  and  OPNAVINST  9072.2  [7]. 

With  the  growth  in  high  performance  computing  power 
and  rising  costs  of  live-fire  testing  for  full  ship  shock 
trials,  the  U.S.  Navy  has  pushed  for  the  development  of 
computer  simulation  techniques  to  supplement  and  enhance 
the  shock  testing  of  new  platforms.  Working  in  conjunction 
with  a  German  Defense  Contractor,  Industrieanlagen- 
Betriegsgesellschaf t  mbH  (IABG),  and  Lawrence-Livermore 
National  Laboratory  the  U.S.  Navy  developed,  fielded,  and 
validated  the  Dynamic  System  Mechanics  Advanced  Simulation 
(DYSMAS)  hydrocode  to  model  the  fully-coupled,  fluid- 
structure  interaction  problem  of  an  UNDEX  event  on  a  ship. 

The  recent  focus  on  naval  operations  in  littoral 
waters,  coupled  with  the  delivery  of  the  Littoral  Combat 
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Ship  (LCS) ,  presents  unique  challenges  to  the  field  of 
shock  testing  and  simulation.  As  nearly  all  UNDEX  testing 
has  been  conducted  in  deep  water,  the  effect  a  shallow 
water  environment  has  on  the  UNDEX  response  of  a  ship  has 
not  been  extensively  investigated.  The  sinking  of  the  Ex- 
Lutjens  in  the  Baltic  Sea  in  is  one  of  the  few  documented 
tests  in  littoral  waters  in  which  DYSMAS  was  used  to 
simulate  the  UNDEX  event.  The  simulation  took  into  account 
the  effect  of  the  ocean  bottom  through  the  application  of 
an  equation  of  state  in  the  fluid  solver  of  DYSMAS,  which 
was  used  to  represent  the  bottom  soil  [8]  .  This  method  of 
bottom  modeling  utilizing  a  fluid  with  representative 
properties  and  behavior  of  soils  has  been  thoroughly  tested 
and  validated  [9],  [10],  [11]. 

B.  SCOPE  OF  RESEARCH 

At  its  most  fundamental  level,  the  accepted  bottom 
modeling  method  used  in  DYSMAS  treats  the  soil  as  an 
inviscid  fluid.  This  research  will  develop  an  alternative 
approach  that  utilizes  the  structural  solver  of  DYSMAS  to 
model  the  bottom  as  a  solid,  finite  element  structure.  This 
Lagrangian  solid  bottom  modeling  approach  will  be  compared 
to  the  current  bottom  modeling  technique  to  determine  its 
validity  and  potential  benefits.  The  Lagrangian  solid 
bottom  modeling  method  provides  the  capability  to  model 
contoured  bottom  profiles,  while  the  current  method  can 
only  create  horizontal  surfaces.  A  parametric  study  of  this 
capability  will  be  conducted  to  determine  the  effect  that 
contoured  bottom  profiles  have  on  the  response  of  a  ship 
subjected  to  an  UNDEX  event  in  littoral  waters. 
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II.  UNDERWATER  EXPLOSION  THEORY 


A.  SHOCKWAVE  OF  INITIAL  EXPLOSION 

The  first  of  two  major  aspects  of  an  underwater 
explosion  (UNDEX)  is  the  high  pressure  "shock"  wave  caused 
by  the  detonation  of  a  high  explosive,  such  as 
trinitrotoluene  (TNT)  .  In  the  case  of  TNT,  the  pressure 
difference  between  the  explosive  and  the  surrounding  water 
is  on  the  order  of  2xl06  psi  (~15  GPa)  .  This  large  pressure 
difference  causes  the  water  to  compress  and  the  wave 
initially  propagates  approximately  five  times  faster  than 
the  acoustic  speed  in  water  (^5,000  ft/sec).  The  pressure 
wave  as  a  function  of  time  passing  through  an  arbitrary 
point  is  steep-fronted,  with  a  backside  that  decays 
exponentially  over  a  few  milliseconds. 

1.  Acoustic  Wave  Assumption 

The  initial  wave  velocity  can  be  five  times  larger 
than  the  acoustic  speed.  Researchers  have  shown  that  the 
wave  quickly  slows  to  acoustic  speed  as  the  wave  front 
pressure  difference  decreases  radially  outward.  Therefore 
the  shockwaves  are  approximated  as  acoustic  plane  waves 
with  small  compressions  and  speeds.  This  assumption  allows 
the  derivation  of  the  Equation  (2.1),  which  relates  added 
fluid  pressure  ( P )  from  the  shockwave  to  the  wave  velocity 
(u)  [2]. 

P  =  P0C0U  (2.1) 

The  mass  density  (pG)  and  the  acoustic  wave  speed  (CG) 
of  the  fluid  are  assumed  to  be  constant  properties.  This 
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relationship  assumes  the  fluid  is  uniform,  compressible, 
and  nonviscous.  The  acoustic  wave  assumption  validation  is 
demonstrated  in  Figure  1 . 


Figure  1.  Shock  Wave  Pressure  Distribution  of  300  lbf 

TNT  Charge  at  Three  Times.  From  [2] 

Figure  1  shows  the  magnitude  of  the  pressure  wave  as 
it  passes  through  three  distances  (R)  from  the  epicenter  of 
the  explosion.  The  solid  lines  represent  the  empirical 
shockwave  values.  The  dashed  lines,  which  lag  the  shock 
wave  by  less  than  five  feet  in  Figures  1  (b)  and  1  (c)  ,  are 
the  acoustic  approximations  taken  at  the  same  time  as  the 
associated  solid  line  shockwave.  With  the  acoustic  speed  in 
water  at  five  feet  per  millisecond,  it  is  thus  safe  to  make 
the  acoustic  wave  assumption.  Cole  has  shown  this 
approximation  to  be  valid  for  radii  from  the  charge  that 
are  between  10  and  100  times  the  radius  of  the  charge. 

2 .  Behavior  at  Interfaces 

Shockwave  behavior  at  an  interface  such  as  the  water- 
seabed  or  the  air-water  is  governed  by  a  combination  of 
Snell's  law  and  the  continuity  of  pressure  and  velocity  at 
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the  interface.  An  arbitrary  fluid  interface  is  diagramed  in 
Figure  2  and  will  be  the  basis  of  the  equations  that 
follow . 


\ 


Figure  2 .  Arbitrary  Interface  Geometry 


Snell'' s  Law  demonstrates  the  relationship 
medium's  acoustic  wave  speed  (C)  and  each  wave's 
incidence  with  reference  to  a  normal  vector 
interface  ( a )  to  one  another. 
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Applying  continuity  of  pressure  and  velocity,  along 
with  Equations  (2.1)  and  (2.2),  to  the  interface  yields  the 
following  relationship  between  the  pressure  and  velocity  of 
the  incident,  reflected,  and  transmitted  shockwaves. 
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(2.4) 


Applying  these  relationships  at  a  rigid  boundary  ( P2C2 
»  P1C1)  yields  a  compressive  wave  reflection  with  the 

magnitudes  of  both  the  incident  and  reflected  waves 
pressure  and  velocity  being  egual .  By  definition  the 
velocity  of  a  rigid  boundary  is  zero. 

Alternatively,  the  shockwave  behavior  at  a  free 
surface,  such  as  the  air-water  interface,  ( P2C2  «  P1C1) 

causes  the  creation  of  a  tensile  reflection  or  rarefaction 
wave.  It  is  referred  to  as  a  tensile  wave  because  Equation 
(2.3)  gives  the  result  that  Prefi  equals  the  negative  of  Pinc- 
The  subsequent  interactions  of  the  incident  (compressive) 
and  rarefaction  (tensile)  waves  can  lead  to  the  phenomenon 
of  bulk  cavitation,  which  will  be  further  discussed  later 
in  this  chapter.  Whether  or  not  bulk  cavitation  occurs,  the 
passage  of  a  compressive  and  tensile  wave  through  a  point 
at  different  times  causes  a  unique  effect  called  the 
pressure  cutoff.  This  effect  is  demonstrated  in  the 
pressure-time  history  in  Figure  3.  The  time  between  the 
passing  of  the  two  waves  is  called  the  cut-off  time.  The 
net  effect  of  the  cut-off  on  a  point  is  an  abrupt  halt  to 
the  impulse  on  a  structure  at  that  point.  Figure  4  provides 
a  visual  depiction  of  the  UNDEX  environment  with  its 
associated  waves. 
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P(t)  of  Initial  Shockwave  for  W=60lbf  of  HBX-1 


Figure  3.  MATLAB  Figure  of  Pressure  History 


Figure  4 .  UNDEX  Geometry 
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3 .  Empirical  Formulas 


Empirical  formulas  to  approximate  several  important 
properties  of  UNDEX  events  have  been  developed  by  Cole  and 
others.  These  equations  were  fit  from  the  interdependence 
of  charge  weight  (i¥)  in  pounds  and  wave  propagation 
distance  (R)  in  feet  for  a  spherical  charge.  The  two  most 
important  formulas  are  for  the  peak  pressure  (Pmax)  in 
pounds  per  square  inch  (psi)  and  the  decay  constant  (6)  in 
milliseconds.  Ki,  Ai,  K2,  and  A2  are  experimental  constants, 
which  depend  on  the  explosive  material  type  [2] . 
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max 


K , 


1  *  J 
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(2.5) 

(2.6) 


B .  GAS  BUBBLE 

After  the  emission  of  the  shockwave,  the  remains  of 
the  exploded  and  unexploded  gas  are  contained  within  a  high 
pressure  gas  bubble  whose  behavior  dominates  the  fluid 
environment  following  the  dissipation  of  the  initial  shock 
wave . 


1 .  Bubble  Motion 

With  the  pressure  of  the  bubble  much  greater  than  the 
surrounding  hydrostatic  pressure,  the  bubble  expands 
rapidly  in  a  roughly  spherical  shape.  The  expansion 
continues  past  the  point  at  which  bubble  pressure  equals 
the  hydrostatic  pressure.  This  is  due  to  the  outward 
momentum  of  the  water.  Once  the  bubble  reaches  its  maximum 
size,  the  surrounding  water  pressure  begins  to  collapse  the 
bubble.  At  a  certain  point,  the  gas  products  inside  the 
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bubble  are  compressed  to  the  point  of  igniting  the 
unexploded  gas.  This  second  explosion  emits  a  shockwave 
known  as  the  first  bubble  pulse.  This  process  of  expansion 
and  contraction  continues  until  the  energy  of  the  bubble  is 
expended  or  the  bubble  "vents"  to  the  atmosphere. 

Buoyancy  and  drag  forces  act  upon  an  oscillating 

bubble  as  it  rises  in  the  water.  The  speed  at  which  the 
bubble  rises  is  inversely  proportional  to  the  bubble's 
diameter.  When  the  bubble  diameter  is  at  a  maximum  the  drag 
surface  area  is  at  its  maximum.  When  the  bubble  contracts 

its  vertical  velocity  increases  until  it  begins  to  expand 

again.  While  it  is  apparent  that  due  to  the  gas'  buoyancy 
the  bubble  must  rise,  the  proximity  to  underwater 

structures  can  significantly  impair  this  movement  [3] . 
Figure  5  details  the  timing  relationship  between  bubble 
size  and  pulses. 


Figure  5.  Bubble  Phenomena  of  UNDEX.  From  [3] 
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2. 


Bubble  Effects 


The  two  major  effects  of  the  bubble  on  ships  are 
follow-on  pulses  and  the  inducement  of  hull  whipping.  Hull 
whipping  occurs  when  the  period  of  the  bubble  pulse 
coincides  with  the  ship's  lowest  natural  frequency.  The 
effect  of  the  pulse  shockwaves  is  significantly  lower  than 
that  of  the  initial  explosion.  Figure  6  is  a  summary  of 
Aron's  energy  partition  of  an  UNDEX  event.  Within  the 
percentage  for  each  shockwave,  roughly  half  is  imparted  as 
the  acoustic  wave.  The  remainder  is  spent  through  flow 
dissipation  and  other  losses  [1]  . 


■  Initial  Shockwave 

■  1st  Pulse 

■  2nd  Pulse 


Figure  6.  Energy  Partition  of  Bubble  Shockwaves 

3 .  Bubble  Formulas  and  Correction  Factors 

While  there  are  various  empirical  methods  of  modeling 
gas  bubble  oscillations,  this  work's  focus  on  shallow 
waters  will  limit  the  scope  bubble  analysis.  In  order  to 
accurately  model  the  extent  and  duration  of  impact  the 
bubble  will  have  on  the  simulations,  an  approximation  of 
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the  bubble's  maximum  radius  and  period  of  the  first  pulse 
is  necessary.  Due  to  the  proximity  of  the  air-water  and 
water-seabed  interfaces  in  shallow  waters,  a  correction 
factor  must  be  taken  into  account. 


Following  the  same  parameters  for  weight  and  distance 
as  Equations  (2.5)  and  (2.6),  the  maximum  bubble  radius 
(Amax)  in  feet  is  [2]  : 


^max  -^6 


w 


D  +  33J 


(2.7) 


Equation  (2.8)  is  used  to  determine  the  bubble  period 
(T)  in  seconds  with  a  correction  factor  for  surface 
proximity.  The  numerical  constant  (a)  is  equal  to  0.1  when 
Amax/D  is  less  than  0.5.  When  Amax/D  is  greater  than  0.5, 
then  a  is  approximated  between  0.1  and  0.2  [12] . 


T 


(2.8) 


C .  CAVITATION 


Cavitation  is  the  creation  and  subsequent  collapse  of 
vapor  bubbles  in  a  fluid,  which  is  the  result  of  the 
pressure  dropping  below  the  fluid's  vapor  pressure.  UNDEX 
events  provide  the  possibility  for  two  distinct  areas  of 
cavitation . 


1 .  Bulk  Cavitation  Zone 

As  discussed  previously,  bulk  cavitation  occurs  when  a 
tensile  wave  lowers  the  total  pressure  at  a  point  below  the 
vapor  pressure  of  the  fluid.  The  result  is  a  random 
distribution  of  small  bubble  formation  throughout  the 
cavitation  region.  The  horizontal  extent  of  the  cavitation 
zone  is  typically  an  order  of  magnitude  larger  than  the 
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depth  to  which  it  extends.  The  bulk  cavitation  zone 
subsequently  collapses  or  "zips  up"  shortly  after  the 
initial  shock  wave  has  diminished.  This  collapse  emits  a 
shockwave  of  its  own,  which  is  much  less  than  the  initial 
shockwave,  but  is  noticeable  in  the  structural  response. 
Figure  7  depicts  the  radial  and  vertical  dimensions  of  the 
bulk  cavitation  zone  for  a  60  pound  charge  of  HBX-1  at  a 
depth  of  32.8  feet.  Equations  (2.9)  and  (2.10),  along  with 
the  variables  listed  in  Table  1,  are  the  empirical  formulas 
used  to  determine  the  upper  and  lower  bounds,  respectively, 
of  the  bulk  cavitation  zone  [13] . 


Figure  7.  Bulk  Cavitation  Zone  for  60  lbf  of  HBX-1 
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Table  1.  Bulk  Cavitation  Variables 


Symbol 

Variable 

Units 

X 

Lateral  distance  of  position  from 
Charge 

ft 

Y 

Depth  of  position  from  free 
surface 

ft 

D 

Depth  of  Charge  from  free  surface 

ft 

ri  = 

1  1 

CN 

K 

+ 

Cl 

Rs 

1 

0.5 

Path  from  Charge  to  position 

ft 

r2  = 

( D  +  y )  +x2 

0.5 

Path  from  Image  Charge  to 
position 

ft 

w 

Weight  of  Charge 

lbf 

c 

Acoustic  Velocity  of  Fluid 

f t/msec 

e 

Decay  Constant  from  Equation 
(2.6) 

msec 

V 

Specific  Weight  of  Fluid 

lbf/  ft3 

Pine 

Pressure  of  Incident  Wave 

lbf/  ft2 

Ki,  A  i,  K2,  A  2 

Coefficients  of  Charge  Material 

N/A 

2 .  Local/Hull  Cavitation 

In  addition  to  bulk  cavitation,  localized  cavitation 

occurs  in  the  area  immediately  surrounding  a  ship's  hull. 
This  effect  is  best  described  in  its  initial  stages  by 

Taylor  Plate  Theory. 

a.  Taylor  Plate  Theory 

This  theory  is  a  1-D  examination  of  a  fluid-plate 

interface  in  order  to  understand  the  structural  response  of 

the  plate  due  to  pressure  wave  impact  through  the  water.  A 

depiction  of  the  interface  is  illustrated  in  Figure  8 . 
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Equation  (2.11)  is  the  equation  of  motion  for  the  plate. 
The  primary  assumption  of  this  derivation  is  that  the 
fluid-structure  interface  (FSI)  is  maintained  at  all  times. 


mw(t)  +  [Jb  +  p0C0\  w(t)  +  kw(t)  =  q(t)  -  2Pinc(0,  t)  (2.11) 
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Figure  8 . 


Taylor  Plate  Geometry.  From  [14] 


Jb.  Separation 

At  some  time  after  the  impact  of  the  pressure 
wave,  depending  on  the  density  of  the  plate,  the  plate  may 
be  moving  at  a  speed  faster  than  the  water  can  maintain. 


These 

differences 

in  speed 

create  a 

tensile 

pull 

on 

the 

fluid 

at  the  FSI. 

At  some 

point  the 

pressure 

may 

drop 

at 

the  FSI  below  vapor  pressure  of  the  fluid  and  local 
cavitation  will  occur. 

Throughout  the  simulations,  as  the  structure  is 

vibrating,  the  presence  of  local  cavitation  in  the  areas 

immediately  around  the  hull  is  expected.  Additionally, 

since  the  intent  of  this  research  is  to  model  the  soil 
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bottom  as  a  structure,  which  can  vibrate  depending  on  the 
material  selected,  local  cavitation  might  be  possible  at 
the  soil-water  interface. 

D.  SHOCKWAVE -BOTTOM  SEDIMENT  INTERACTION 

As  discussed  in  Section  II. A. 2,  a  pressure  wave 
incident  on  a  rigid  interface  will  be  reflected  as  a 
compressive  wave  of  equal  pressure  and  speed.  If  the  ocean 
bottom  were  assumed  to  be  a  rigid  body,  it  would  simplify 
the  modeling  requirements.  However,  the  ocean  bottom  can  be 
composed  of  a  range  of  materials  yielding  multiple  shock 
responses.  For  example,  highly-porous ,  bottom  sediment  can 
reflect  a  tensile  (or  rarefaction)  wave  instead  of  a 
compressive  wave  [10] . 

The  two  aspects  of  waves  in  soils  that  have  been  well 
researched  and  tested  are  the  response  of  dry,  porous  soil 
undergoing  air-blast  events  and  the  properties  of 
submerged,  porous  soil  with  respect  to  their  reflection, 
transmission,  and  attenuation  properties  when  subjected  to 
acoustic  waves  at  pressure  levels  significantly  lower  than 
UNDEX  shockwaves.  As  a  result,  most  available  soil  material 
models  used  are  focused  on  these  two  aspects.  The 
particular  shortcoming  is  the  lack  of  bottom  sediment 
experimental  data,  which  is  fitted  to  already  established 
material  models. 

1 .  Soil  Properties 

All  naturally  occurring  soil  is  a  three-phase  system 
consisting  of  solid  particles,  water,  and  air.  The  solid 
particles  form  structures  of  varying  degrees  of  order.  The 
voids  within  this  structure  are  filled  either  by  air  or 


15 


water.  Additionally,  there  is  a  distinct  difference  in 
behavior  between  cohesive  (clay)  and  non-cohesive  (sand) 
soils  [15]  .  In  general,  the  tensile  strength  of  all  soils 
is  orders  of  magnitude  smaller  than  it's  compressive  and 
shear  strengths  [16]. 

2.  In-Situ  Stresses  of  Saturated  Soils 

The  definition  of  a  saturated  soil  is  one  in  which 
water  has  filled  all  of  the  voids  in  the  solid  structure. 
An  understanding  of  the  static  stresses  in  the  soil  is 
necessary  in  order  to  effectively  model  boundary  conditions 
on  the  bottom  model.  The  total  static  stress  at  any  point 
in  the  soil  column  can  be  found  by  applying  Equation 
(2.12)  .  The  equation  is  based  upon  the  assumption  that  the 
water  in  the  soil  is  not  flowing  vertically  through  the 
soil.  Equation  (2.12)  determines  the  total  stress  (a)  at  a 
depth  below  the  air-water  interface  using  the  unit  weights 
of  water  (yw)  and  saturated  soil  (ysat)  r  the  depth  from  the 
air-water  interface  to  the  water-soil  interface  (H)  ,  and 
the  depth  from  the  water-soil  interface  to  the  position  of 
interest  (HA)  [15]  . 


a  =  Hlw  +  (HA  -  H)  7sat 


(2 . 12) 
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III.  MODELING  AND  SIMULATION  USING  DYSMAS 


A.  SOFTWARE  SUITE 

The  software  package  that  was  used  for  this  research 
is  the  Dynamic  System  Mechanics  Advanced  Simulation 
(DYSMAS)  code.  DYSMAS  is  a  fully  coupled,  six  degree-of- 
freedom,  hydrocode  that  is  designed  to  simulate  three- 
dimensional,  UNDEX  events.  The  software  consists  of  three 
major  components.  The  Gemini  software  is  an  Eulerian  solver 
for  the  fluid  environment.  The  structural  solver, 
Dyna_N(3D),  evaluates  the  structural  response  using  a 
Lagrange  method.  While  several  other  structural  solvers  can 
be  used  with  DYSMAS,  only  Dyna_N(3D)  was  used  for  this 
research.  The  final  and  most  important  component  is  the 
Standard  Coupler  Interface  (SCI),  which  passes  information 
between  the  two  solvers  at  the  end  of  each  Eulerian  time 
step  to  maintain  the  fully  coupled  interface. 

B .  GEMINI 

The  Gemini  code  is  a  finite-difference,  Euler  equation 
solver.  It  was  specifically  designed  to  simulate  the  UNDEX 
phenomena  of  shockwaves  and  bubble  jetting.  It  is  capable 
of  solving  flow  fields  involving  several  fluid  types. 

1 .  Theory 

Gemini  solves  the  fluid  mesh  by  running  a  higher  order 
Godunov  method  algorithm  to  solve  the  fluid  domain  at  each 
time  step.  This  algorithm  is  applied  to  each  Euler  cell  in 
a  one-dimensional  approach  through  each  principal  direction 
at  a  time  solving  the  Euler  equations  for  conservation  of 
mass,  momentum,  and  energy.  A  major  assumption  of  the  Euler 
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equations  is  that  the  flow  is  frictionless  or  inviscid. 
This  assumption  of  an  inviscid  fluid  limits  Gemini's 
capability  to  fully  model  a  Lagrangian  solid  bottom  as  most 
soils  can  support  shear  loading  [17]  . 

2 .  Components 

There  are  several  additional  components  to  the  Gemini 
code.  GemGrid,  PreGemini,  Float,  and  Prestress  are 
preprocesses.  GemGrid  allows  the  user  to  setup  up  the 
three-dimensional,  Euler  cell  grid.  This  allows  for  grid 
refinement  around  specific  areas  of  interest  in  the  flow 
field,  specifically  the  fluid  bubble  and  the  ship  hull. 
PreGemini  fills  the  Euler  cell  grid  with  the  user-defined 
material  equations  of  state  and  their  initial  properties 
(energy,  density,  etc.).  The  Float  and  Prestress  processes 
allow  the  user  to  accurately  position  and  prestress  the 
Dyna_N  structure  in  the  fluid  domain  prior  to  the  transient 
analysis.  Following  the  completion  of  the  simulation  by  the 
Gemini  solver,  the  postprocessing  of  the  data  is  completed 
by  GemHis  and  GemField.  GemHis  processes  the  recorded  data 
at  specific  points  of  interest  throughout  the  fluid  for 
analysis.  GemField  generates  contour  plots  of  the  flow 
field  at  specified  times  for  the  fluid  domain  and  interface 
segments.  The  relationship  of  each  Gemini  component  to  the 
overall  code  is  diagramed  in  Figure  9.  The  light  blue  boxes 
indicate  the  various  user-defined  input  files  that  are 
required  [18]. 
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Figure  9. 


Gemini  Data  Flow  Chart.  From  [18] 


C.  DYNA_N(3D) 

Dyna_N  is  a  nonlinear,  explicit,  three-dimensional 
finite  element  code  that  analyzes  dynamic  structural 
responses . 


1 .  Theory 

The  software  solves  a  discrete  formulation  of  the 
linear  second  order  differential  equation  of  motion: 

mx  +  cx  +  kx  =  p{t )  (3.1) 

It  can  also  accommodate  non-linearity  in  the  structure 
with  the  following  equation: 

mx  +  cx  +  fint  =  p  (t)  (3.2) 

Dyna_N  solves  the  combination  of  these  two  equations, 
which  has  a  form  of: 

[M]{x}  +  {F}-{/>}  =  0  (3.3) 
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where  {P}  represents  the  body  forces  and  external  loads, 
{x}  is  the  displacement  vector,  [x''}  is  the  acceleration 
vector,  [M]  is  the  structure's  mass  matrix,  and  {F}  handles 
the  internal  non-linearity  [19]. 

The  solver  then  utilizes  an  explicit  central 
differencing  scheme  to  step  through  each  time.  The  use  of 
the  explicit  time  step  integration  dictates  that  the 
solution  is  only  stable  when  the  time  step  is  kept  below  a 
certain  value,  commonly  known  as  the  Courant  number. 
Essentially,  this  means  that  the  time  step  must  be  kept 
small  enough  to  allow  a  wave  to  propagate  through  the  space 
discretization.  For  example,  if  the  wave  speed  is  5000 
ft/sec  and  the  discretized  length  is  only  one  foot,  the 
time  it  takes  the  wave  to  pass  through  the  element  is  0.2 
milliseconds.  If  the  time  step  were  any  larger,  the  wave 
would  pass  through  the  element  without  the  solver  having 
the  opportunity  to  account  for  it.  In  the  end,  this  means 
that  the  maximum  time  step  is  driven  by  the  smallest 
element  in  the  structure.  The  standard  coupler  interface 
will  choose  the  smaller  of  the  two  time  steps  between 
Gemini  and  Dyna_N. 


2.  Pre-/Post-Processing 

The  preprocessing  of  all  finite  element  models  was 
done  using  the  DYSMAS  Pre-Processor  2010.  The  software 
allows  for  the  creation  of  the  model's  structure, 
assignment  of  material  properties,  boundary  conditions, 
body  forces,  motions,  and  fluid-structure  interface  segment 
definitions.  The  structural  model  is  then  written  into  the 
specific  input  cards  required  Dyna_N  [20] . 
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The  use  of  the  preprocessor  in  generating  the  ocean 
bottom  structure  proved  problematic  as  it  did  not  have  the 
ability  to  define  several  different  material  models, 
boundary  conditions,  and  functions  that  were  tested.  These 
features  were  present  in  Dyna_N,  but  no  user-interface  in 
the  preprocessor  had  been  created.  As  a  result,  several 
features  of  each  model  had  to  be  tediously  entered  into  the 
Dyna_N  input  cards  after  the  preprocessor  had  been  run. 

The  DYSMAS  postprocessor  allows  the  user  to  visualize 
the  entire  fluid-structure  dynamic  response.  It  contained 
features  that  allowed  for  the  extraction  and  analysis  of  a 
wide  range  of  response  data  from  both  the  fluid  domain  and 
the  structure. 

D.  STANDARD  COUPLER  INTERFACE  (SCI) 

The  SCI  is  the  key  component  in  the  DYSMAS  software 
suite  that  allows  the  simulations  to  be  fully  coupled.  The 
SCI  is  required  because  in  order  to  pass  information 

between  the  fixed  Euler  grid  and  a  finite  element  model 

whose  structural  nodes  and  interfaces  move.  Figure  10 
depicts  this  link  as  well  as  the  type  of  information 

passed.  The  main  functions  the  SCI  accomplishes  are  [17] : 

•  Build  a  grid  representation  of  the  structural 
interface  into  the  Euler  mesh 

•  Enforce  boundary  conditions  at  the  FSI 

•  Activate/De-activate  cells  as  the  FSI  passes 

through  the  Euler  mesh 

•  Calculate  nodal  loads. 
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Figure  10.  SCI  Block  Diagram.  From  [18] 

E.  BASIC  DYSMAS  SIMULATION  SETUP 

The  creation  of  each  DYSMAS  simulation  used  in  this 
research  involved  nearly  all  the  same  basic  steps.  These 
included  generating  the  Euler  grid,  filling  the  fluid 
domain  with  the  correct  hydrostatic  conditions,  and 
developing  the  finite  element  model.  For  each  of  these 
steps  several  different  sources  of  best  practices  were 
utilized  as  decision-making  tools. 

1 .  Euler  Grid  Generation 

When  defining  the  Euler  grid,  the  most  important 
determination  was  the  grid  refinement  immediately 
surrounding  the  gas  bubble  and  the  ship  structure. 
According  to  Prendergast,  DYSMAS  simulations  are  highly 
sensitive  to  the  mesh  size  surrounding  the  explosion.  The 

smaller  the  mesh  size,  the  closer 
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in  accuracy  the 


simulation  is  to  both  empirical  predictions  and 
experimental  data  in  maximum  pressure,  bubble  radius,  and 
bubble  period.  Additional  simulations  demonstrated  that 
Gemini's  free  boundary  conditions  did  not  adequately  allow 
fluid  to  flow  out  of  and  back  into  the  fluid  domain  quickly 
enough  if  the  fluid  domain  was  too  small  in  comparison  to 
the  maximum  bubble  diameter.  This  finding  was  particularly 
applicable  to  the  shallow  water  simulations  where  the 
overall  fluid  depth  was  restricted.  Based  on  Prendergast ' s 
findings  and  discussions  of  best  practices  with  Didoszak, 
the  parameters  in  Table  2  were  used  to  determine  the  size 
of  both  the  fine  mesh  and  also  overall  extent  of  the  fluid 
domain  [21],  [ 22 ]  . 


Table  2.  GemGrid  Decision  Matrix 


2-D  Mesh  Size 

2  cm  (x,  z ) 

3-D  Fine  Mesh  Size 

<  Minimum  length  of  smallest 

interface  element  (Default  =  10cm) 

Extent  of  Bubble  Fine 

>  1.71  x  Max.  Bubble  Radius  in 

Mesh 

(x,  y,  z) 

Extent  of  Ship  Fine 

Mesh 

>=  100cm  on  each  side 

Distance  from  Charge 

to  Fluid  Boundary 

=  10  x  Max.  Bubble  Radius  in  (x,y) 

Thickness  of  Air  Layer 

=  1,000cm 
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2. 


Fluid  Domain  Setup 


The  fluid  domain  was  initialized  through  the  Pregemini 
input  file.  The  critical  step  in  this  procedure  was  the 
seeding  of  the  fluid  domain  with  the  correct  materials  at 
the  appropriate  states.  In  order  to  simulate  the  different 
materials  such  as  explosives,  air,  water,  etc.  Gemini 
utilizes  various  Equations  of  State  (EOS) .  These  EOSs  are 
used  to  relate  the  pressure,  energy,  and  density  through 
the  solution  of  the  fluid  environment.  The  EOSs  used  in 
this  research  are  provided  in  Tables  3  through  7  [18] . 


Table  3.  Gemini  EOS  for  Air 


Material : 

Air 

EOS: 

Gamma-Law 

Form : 

P  =  (7  -  1)  Pe 

Properties 

Cavitation  Pressure 

Pcav 

0  . 

dynes /cm2 

Reference  Density 

Po 

0.0013 

g/cm3 

Reference  Specific  Energy 

eD 

1 . 92307  69e+9 

erg/g 

Speed  of  Sound  in 

Cavitation  Region 

Ccav 

3 . 2  8e+4 

cm/s 

Minimum  Density 

P floor 

1  .D-06 

g/cm3 

Minimum  Specific  Energy 

&  floor 

1  .D-4 

erg/g 

Gamma 

Y 

1.4 

N/A 

Table  4.  Gemini  EOS  for  HBX-1  (Unburned) 


Material : 

HBX-1  Unburned 

EOS: 

Tait  for  Unburned  Explosive 

Form : 

P  =  Po  +  B 

7  N 

p  i 

\Po) 

/ 

Properties 

Cavitation  Pressure 

Pcav 

0  . 

dynes /cm2 

Reference  Density 

Po 

1.720 

g/cm3 

Reference  Specific  Energy 

eD 

6 . 520E+10 

erg/g 

Speed  of  Sound  in 
Cavitation  Region 

Ccav 

4 . 0  6et5 

cm/s 

Minimum  Density 

P floor 

1  .E-06 

g/cm3 

Minimum  Specific  Energy 

€ floor 

1  .E-4 

erg/g 

Gamma 

Y 

1 

N/A 

Ambient  Pressure 

Po 

1  .E+6 

dynes /cm2 

Unburned  Explosive 

Density 

Po 

1.72 

g/cm3 

CJ  Density 

Pcj 

2.25567 

g/cm3 

CJ  Specific  Energy 

eCJ 

7 . 455692E+10 

erg/g 

CJ  Pressure 

Pcj 

1 . 355396E+11 

dynes /cm2 

Detonation  Velocity 

D 

5 . 75045E+05 

cm/s 

Factor 

F 

0 . 95 

N/A 

Table  5.  Gemini  EOS  for  HBX-1  (Exploded) 


Material : 

HBX-1  (Exploded) 

EOS: 

JWL 

Form : 

{  uop  )  _Ritr  lop  )  _R2if 

p  =  A  1  —  -  e  p  +  B  1  —  -  e  p  +  c ope 

R,  p„  R?P-, 

lr  o  J  2r  o  / 

Properties 

Cavitation  Pressure 

Pcav 

o . 

dynes /cm2 

Reference  Density 

Po 

1.720 

g/cm3 

Reference  Specific  Energy 

6 . 520Etl0 

erg/g 

Speed  of  Sound  in 

Cavitation  Region 

Ccav 

4 . 0  6e+5 

cm/s 

Minimum  Density 

P floor 

1  .E-06 

g/cm3 

Minimum  Specific  Energy 

^ floor 

1  .E-4 

erg/g 

constant 

A 

5 . 183E+12 

N/A 

constant 

B 

4.39E+9 

N/A 

constant 

Ri 

5.183 

N/A 

constant 

r2 

3.5E-1 

N/A 

y  -  1 

CO 

0.25 

N/A 
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Table  6. 


Gemini  EOS  for  Water 


Material : 

Water 

EOS: 

Tillotson 

Form : 

p  =  Po  +  up  (e  -  e0)  +  A 

(  \ 

-l-i 

,Po 

+  B 

f  \ 

-l-i 

k  1 

2 

+  c 

(  \ 

-l-i 

k  1 

3 

Properties 

Cavitation  Pressure 

Pcav 

50000  . 

dynes /cm2 

Reference  Density 

Po 

1.0 

g/cm3 

Reference  Specific  Energy 

eD 

3 . 5420001E+9 

erg/g 

Speed  of  Sound  in 

Cavitation  Region 

Ccav 

147600.0 

cm/s 

Minimum  Density 

P floor 

9 . 999998E-03 

g/cm3 

Minimum  Specific  Energy 

^ floor 

-9 . 99998E+10 

erg/g 

y  -  1 

CO 

.28000000119 

N/A 

Ambient  Pressure 

Po 

1  .E+6 

dynes /cm2 

constant 

A 

2 . 200000E+10 

N/A 

constant 

B 

9 . 540000E+10 

N/A 

constant 

C 

1 . 457E+11 

N/A 

Table  7 .  Gemini  EOS  for  Clay 


Material : 

Clay 

EOS: 

Mie-Gruneisen 

Form : 

P  =  roPo 

P 

with : 

e  - 

=  1 

Po  Cs2 

rwvl 

-7  r  = 

"  and  k=0 

d 

+  P°  — 

dp 

(  7  ) 

cs2y 

1(1  -  5m)  1 

+  Po 

Properties 

Cavitation  Pressure 

Pcav 

20000  . 

dynes /cm2 

Reference  Density 

Po 

2.023 

g/cm3 

Reference  Specific  Energy 

e0 

0  . 

erg/g 

Speed  of  Sound  in 

Cavitation  Region 

Ccav 

260000.0 

cm/s 

Minimum  Density 

P floor 

1.0E-02 

g/cm3 

Minimum  Specific  Energy 

^ floor 

-9 . 99998E+99 

erg/g 

Gamma 

rQ 

0 . 97 

N/A 

Slope  of  Shock  Speed  - 
Particle  Velocity  Plot 

s 

1 .86 

N/A 

Square  of  Reference  Sound 
Speed 

c  2 

4 . 048E+10 

cm2/s2 
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3 .  Charge  Parameters 

In  all  simulations,  except  where  specifically  noted,  a 
standard  charge  of  HBX-1  was  used.  The  charge  weight  was  60 
pounds  or  27.215  kilograms.  It  was  placed  at  a  depth  of 
32.8  feet  or  ten  meters.  As  noted  in  Section  II.  A.  3,  by 
keeping  the  charge  weight  and  depth  constant,  the  resulting 
shockwave  and  bubble  parameters  remain  constant.  This 
allows  the  effects  of  the  presence  and  properties  of  the 
Lagrangian  solid  bottom  to  be  determined. 

4 .  Structural  Finite  Element  Model 

The  Floating  Shock  Platform  (FSP)  finite  element  model 
was  used  for  all  simulations  involving  a  floating 
structure.  The  FSP  is  a  well-defined  model  whose  properties 
and  behaviors  under  shock  loading  have  been  thoroughly 
tested  and  validated.  Additionally,  its  relatively  small 
number  of  elements  (12,792)  and  total  size  (approximately 
9m  x  5m  x  2.5m)  allowed  for  quick  simulation  times.  The 
particular  FSP  model  used  has  been  slightly  modified  from 
the  true  structure.  A  thin,  lightweight,  highly  flexible, 
elastic  roof  has  been  attached  over  the  open-topped  FSP. 
This  was  done  in  order  to  be  able  to  define  the  entire  FSP 
as  a  one  singly-wetted  interface  on  all  six  sides. 

In  order  to  accurately  model  the  Lagrangian  solid 

bottom,  the  selection  of  the  appropriate  constitutive 

equation  for  the  material  is  important.  The  most  applicable 

soil  models  available  in  Dyna_N  are  material  types  16,  45, 

and  65.  These  models  were  created  to  model  concrete  and 

geologic  materials.  Due  to  the  complexities  of  these  soil 

models,  in  conjunction  with  the  creation  of  a  new  method  of 

Lagrangian  solid  bottom  modeling,  these  material  types  were 
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not  considered  in  this  research.  In  order  to  ensure  the 
fewest  errors  while  testing  new  modeling  methods,  an 
elastic  material  model  was  utilized.  The  properties  of  the 
material  used  were  developed  from  Bangash' s  research  into 
explosion  dynamics  of  numerous  soil  materials  and  are 
listed  in  Table  8  [23] . 

Table  8.  Lagrangian  Solid  Bottom  Elastic  Soil  Properties 


Property 

Value 

Units 

Density 

2 . 0 

g/  cm3 

Elastic  Modulus 

2  .  OE+11 

dyne/ cm2 

Poisson''  s  Ratio 

0.3 

N/A 

5 .  Summary  of  Simulation  Setup 

All  of  the  values,  properties,  and  decisions  made  in 
the  setup  for  each  simulation  started  from  the  guidance  set 
forth  in  this  section.  This  basic  construct  was  used  to 
develop  and  test  the  initial  Lagrangian  solid  bottom  model. 
Various  modifications  were  made  as  the  model  and  its 
parameters  were  proven.  Figure  11  provides  a  visualization 
of  the  fluid  domain  along  with  the  FSP  used  in  this 
standardized  setup.  A  complete  index  of  simulations  that 
were  conducted  in  this  research  can  be  found  in  the 
Appendix . 
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IV.  LAGRANGIAN  SOLID  BOTTOM  MODELING  IN  DYSMAS 


There  are  two  fundamentally  different  approaches  to 
modeling  a  Lagrangian  solid  bottom  in  DYSMAS.  The  first 
approach  is  depicted  in  Figure  12.  This  method  places  the 
entire  structure  within  the  bounds  of  the  fluid  domain. 
This  approach  closely  resembles  the  placement  of  a  ship' s 
FEM  in  the  fluid  domain.  As  such,  it  would  be  relatively 
easy  to  implement.  However,  the  constitutive  equations  that 
govern  the  behavior  of  the  soil  material  in  both  Dyna_N  and 
Gemini  can  never  be  perfectly  matched. 


Wall 


Figure  12.  Bottom  Fully  Contained  within  Fluid 

The  second  approach  is  to  place  the  Lagrangian  solid 
bottom  only  partially  within  the  fluid  domain  as  shown  in 
Figure  13.  This  allows  the  soil  to  be  treated  as  a  semi¬ 
infinite  domain  and  reduces  the  amount  of  cells  in  the 
fluid  domain.  The  reduction  in  number  of  fluid  cells  in  the 
Euler  grid  has  the  potential  to  decrease  the  computational 
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time  and  resources  needed  for  each  simulation.  This 
approach  requires  the  investigation  of  the  parameters 
involved  in  implementing  the  first  approach  plus  additional 
parameters  unique  only  to  this  method.  For  this  reason,  the 
second  approach  was  chosen  as  the  primary  method  to 
Lagrangian  solid  bottom  modeling  in  this  research. 


NRB 


Figure  13.  Bottom  Outside  the  Fluid 

A.  SIZE  AND  POSITION  OF  LAGRANGIAN  SOLID  BOTTOM 

The  first  consideration  to  be  made  when  developing  the 
ocean  bottom  model  for  Dyna_N  is  the  size  of  the  model. 
Both  the  overall  size  and  individual  element  size  have 
significant  consequences  for  the  simulation's  outcome. 

1 .  Lateral  Dimensions  and  Position 

The  overall  size  is  governed  by  two  factors.  The  first 
is  the  coupling  interface  requirement  of  the  fluid  code, 
Gemini.  Gemini  allows  nodes  describing  surface  elements  and 
their  associated  singly-wetted  interfaces  (SWI)  to  be 
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positioned  on  the  boundary  of  the  fluid  domain.  A  SWI 
occurs  when  the  element  is  only  in  contact  with  the  fluid 
on  one  of  its  sides.  The  edges  of  the  SWI  cannot  terminate 
within  the  fluid  domain.  They  must  entirely  enclose  the 
structure,  or  they  must  extend  to  or  beyond  the  boundaries 
of  the  fluid  domain.  The  degree  of  nodal  position  precision 
required  by  the  interface  requirement  is  within  10-6 
centimeters.  Most  structural  finite  element  pre-processors 
do  not  achieve  that  degree  of  precision.  The  result  is  that 
the  Lagrangian  solid  bottom  must  be  expanded  in  the  lateral 
dimensions  (X,Y)  to  a  size  just  greater  than  the  lateral 
dimensions  of  the  fluid  domain.  In  doing  so,  the  only 
interface  required  is  the  top  surface  of  the  Lagrangian 
solid  bottom.  The  size  requirement  places  this  approach  at 
a  disadvantage  when  compared  with  a  Lagrangian  solid 
bottom,  which  is  wholly  contained  within  the  fluid,  as 
shown  in  Figure  12  [18] . 

The  portion  of  the  initial  shockwave  that  is  reflected 
toward  the  target  by  the  bottom  only  impacts  the  bottom  in 
a  finite  area  between  the  charge  and  the  target.  Leveraging 
this  knowledge,  a  Lagrangian  solid  bottom  entirely  within 
the  fluid  could  be  sized  to  cover  only  the  area  that  would 
reflect  the  shockwave.  This  approach  would  allow  the 
Lagrangian  solid  bottom  to  be  considerably  smaller  than  the 
Lagrangian  solid  bottom,  which  extends  past  the  lateral 
dimensions  of  the  fluid  domain.  The  reduction  in  resources 
required  to  solve  the  Lagrangian  solid  bottom  response  is 
significant.  Given  the  example  simulation  parameters  listed 
in  Table  9,  modeling  the  Lagrangian  solid  bottom  only 
within  the  region  of  interest  would  reduce  the  number  of 
Lagrange  elements  by  99.5%.  This  second  approach  would 
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require  the  addition  of  fluid  cells  to  make  up  for  the 

99.5%  decrease  in  volume.  Thus  the  approach  simply  reduced 
the  computational  requirements  of  the  structural  solver  at 
the  expense  of  increased  computational  requirements  of  the 
fluid  solver.  If  the  bottom  were  made  up  of  two  different 
materials,  an  Euler  and  a  Laqranqe  soil,  the  bottom 

reflection  cannot  be  assumed  to  be  accurate  across  the 

entire  fluid  domain.  This  is  due  to  the  fact  that  the 

equations  of  states  for  the  fluid  and  the  constitutive 
equations  for  the  structure  are  not  perfectly  matched. 
Therefore,  modelinq  the  bottom  only  in  the  reqion  of 
interest  is  unlikely  to  give  a  valid  response  when  compared 
to  empirical  data. 


Table  9.  Bottom  Modeling  Approach  Comparison 


Charge 

Fluid 

Tarq 

ret 

Type 

HBX-1 

X  Width 

84  m 

Position 

Surface 

Weight 

60  lbf 

Y  Width 

84  m 

Separation 

6  m 

Depth 

10  m 

Z  Depth 

35  m 

Bottom  Modeling  Approach 

Large 

Focused 

Shape 

Parallelepiped 

Shape 

Parallelepiped 

X  Width 

84.02  m 

X  Width 

6  m 

Y  Width 

84.02  m 

Y  Width 

6  m 

Thickness 

5  m 

Thickness 

5  m 

#  Elements 

35,280 

#  Elements 

180 

2 .  Vertical  Dimensions  and  Position 

The  vertical  position  of  the  Lagrangian  solid  bottom 
is  obviously  determined  by  the  required  bottom  depth  of  the 
simulation.  However,  the  top  surface  of  the  bottom  cannot 
be  simply  joined  to  the  bottom  of  the  fluid  domain.  In 
Figure  13,  with  the  Lagrangian  solid  bottom  outside  of  the 

fluid,  the  slight  overlap  of  the  fluid  domain  and 
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Lagrangian  solid  bottom  is  required.  The  bottom  is  expected 
to  deform  under  hydrostatic  pressure,  shockwave  pressure, 
and  its  own  self-weight.  The  deformation  is  assumed  to  be 
in  the  negative  Z-direction.  Remembering  the  SWI 
requirement  that  the  interface  cannot  be  more  than  10”6  cm 
from  the  edge  of  the  fluid  domain,  the  requirement  to 
overlap  the  fluid  domain  and  the  bottom  becomes  apparent. 
This  overlap  prevents  the  interface  from  leaving  the  fluid 
domain  due  to  the  deformation.  The  necessity  to  include 
gravity  in  the  structural  solver  requires  the  application 
of  vertical  constraints  on  the  bottom  surface  of  the 
Lagrangian  solid  bottom.  These  constraints  keep  the 
structure  from  bodily  translating  out  of  the  domain.  In  all 
simulations  in  this  research  the  minimum  overlap  was  one 
meter  at  the  deepest  point  of  the  water-soil  interface. 

An  understanding  of  the  required  thickness  of  the 
Lagrangian  solid  bottom  was  developed  through  the 
examination  of  wave  propagation  in  the  solid.  When  a 
shockwave  is  incident  on  an  interface  between  dissimilar 
materials  the  wave  behaves  under  Snell's  law.  The  result  is 
the  possibility  of  a  reflection  and/or  a  transmission  wave. 
The  reflected  wave  is  the  sent  back  into  the  original 
material,  while  the  transmitted  wave  propagates  through  the 
new  material.  When  examining  shallow  water  UNDEX,  the 
properties  and  effect  of  the  initial  reflected  wave  are  the 
primary  interest.  Conversely,  the  wave  transmitted  into  the 
soil  is  of  little  consequence  because  it  is  assumed  to  have 
propagated  through  the  seabed  and  dissipated  over  some 
distance.  Due  to  the  finite  nature  of  the  bottom  model, 
this  effect  is  not  seen.  As  the  transmitted  wave  propagates 

through  the  solid  material,  the  nodal  constraints  that 

35 


prevent  structural  translation  reflect  this  wave  back  up 
through  the  soil  structure.  Upon  reaching  the  water-soil 
interface,  a  shockwave  is  transmitted  into  the  fluid  domain 
that  follows  close  behind  the  initial  bottom  reflection. 
For  the  extent  of  this  research,  this  artificial  wave  is 
referred  to  as  the  retransmission  wave.  The  time  delay 

between  the  bottom  reflection  and  the  retransmission  is 
governed  by  the  thickness  of  the  soil  structure.  As  the 
thickness  increases  so  too  does  the  delay.  In  order  to 
decrease  the  effect  of  the  retransmission,  the  thickness  of 
the  Lagrangian  solid  bottom  must  be  increased  to  a 
relatively  large  size.  The  obvious  limit  to  expanding  the 
thickness  is  once  again  the  limit  of  computational 
resources.  Further  discussion  of  the  retransmission  effect 
and  methods  to  diminish  it  are  discussed  in  greater  detail 
in  Section  IV. B. 2. 

3 .  Element  Size 

There  are  no  set  restrictions  on  the  size  of 

individual  solid  elements  within  the  Lagrangian  solid 
bottom.  All  of  the  following  suggested  guidelines  were 

formed  through  discussions  with  various  subject  matter 
experts  and  review  of  cautions  located  in  the  Dyna_N  User' s 
Manual . 

The  U.S.  Bureau  of  Reclamation  has  done  extensive 

simulations  of  UNDEX  shock  loading  on  dam  structures 
utilizing  the  DYSMAS  software  package.  The  dam  structure's 
finite  element  models  were  constructed  mainly  with  solid 
elements  with  the  material  properties  of  concrete.  An 
examination  of  the  size  of  their  solid  elements  provided  a 
guide  for  the  dimensions  that  were  adopted  in  this 
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research.  In  the  region  of  interest  on  the  dam  structure, 
the  solid  elements  were  approximately  two  to  three  feet  on 
each  side.  The  mesh  expanded  toward  the  bounds  of  the  model 
with  some  elements  reaching  20  to  30  feet  on  a  side.  This 
research  adopted  solid  elements  that  were  one  meter  in  each 
lateral  dimension.  This  lateral  sizing  was  maintained  for 
every  Lagrangian  solid  bottom  model  used  [24] . 

Several  of  the  simulations  used  Lagrangian  solid 
bottoms  in  which  the  profile  changed  vertically  across  the 
fluid  domain.  The  vertical  dimensions  of  the  solid  elements 
were  allowed  to  expand  and  contract  as  necessary  in  order 
to  create  these  contours.  The  only  constraint  placed  upon 
the  vertical  dimensions  was  that  it  must  be  no  more  than 
three  times  as  large  or  as  small  as  the  lateral  dimensions. 
This  allowed  the  vertical  dimension  to  vary  between  three 
meters  and  one  third  of  a  meter  [25] . 

B.  BOUNDARY  CONDITIONS  AND  CONSTRAINTS 

Typical  DYSMAS  simulations  of  UNDEX  events  involve  a 
buoyant  ship' s  finite  element  model  that  is  wholly 
contained  within  the  fluid  domain.  As  discussed  in  Section 
III.D,  the  linkage  between  the  two  is  defined  through  the 
use  of  the  interfaces  and  the  SCI  software.  Through  these 
interfaces  the  shock  loading  and  the  buoyant  forces  are 
applied  to  the  structure.  The  buoyant  forces  on  the 
structure  serve  as  the  boundary  conditions  for  the 
structure . 

Unlike  the  ship,  the  non-buoyant  Lagrangian  solid 

bottom  requires  the  application  of  boundary  conditions 

other  than  the  Gemini  interfaces.  This  requirement  holds 

true  no  matter  whether  the  model  is  wholly  within  the  fluid 
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domain  or  not.  If  there  were  no  boundary  conditions  on  the 
Lagrangian  solid  bottom,  the  gravity  loading  and 
hydrostatic  pressure  loading  would  cause  the  bottom  to 
"fall  away"  from  the  bottom  of  the  fluid  domain. 

1 .  Fixed  Nodal  Displacements 

The  initial  solution  used  was  the  application  of 
Dirichlet,  or  fixed,  boundary  conditions  to  the  model.  All 
the  nodes  on  the  lateral  and  bottom  faces  of  the  structure 
were  constrained  in  all  six  nodal  degrees  of  freedom: 
translation  in  and  rotation  about  the  X,  Y,  and  Z  axes. 
Table  10  details  the  pertinent  simulation  and  model  data 
used  to  simulate  the  effect  of  the  nodal  constraints.  These 
initial  investigations  all  used  a  purely  elastic  bottom 
material  to  avoid  any  material-specific  errors. 


Table  10.  Input  Parameters  for  Run  ID  4-01 


Charge 

Fluic 

1 

Target 

Type 

HBX-1 

X  Width 

84  m 

Position 

N/A 

Weight 

60  lbf 

Y  Width 

84  m 

Separation 

N/A 

Depth 

10  m 

Z  Depth 

35  m 

Bottom 

Dimensions 

Material 

Shape 

Parallelepiped 

Type 

Elastic 

X  Width 

84.02  m 

Density 

2 . 0  gram/ cm3 

Y  Width 

84.02  m 

Elastic 

Modulus 

2  x  1011 
dyne/ cm2 

Thickness 

5  m 

Poisson' s 
Ratio 

0.3 

#  Elements 

35,280 

EOS  Type 

None 

Boundary  Conditions 

All  nodes  on  lateral  faces  and  underside  of  the 
Lagrangian  solid  bottom  were  constrained  from 
translation  and  rotation  in  all  three  axes. 
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This  initial  simulation  using  a  Lagrangian  solid 
bottom  with  nodal  constraints  was  compared  to  a  fluid-only 
simulation  of  the  charge  and  fluid  parameters  that  utilized 
a  perfectly  reflective  bottom  boundary.  Figure  14  compares 
the  pressure  distributions  through  the  fluid  domain  of  each 
simulation  at  30  milliseconds.  Both  simulations  show  a 
similar  bottom  reflection  of  the  initial  shockwave.  The 
difference  is  the  presence  of  a  second  pressure  wave 
propagating  from  the  bottom  in  Figure  14  (B)  . 


RunID  1-00:  Bottom  Wall  Boundary  Condition 
Euler:  Pressure  (Pa) 

Lagrange:  Vertical  Stress  (dyne/cmA2) 

Time  :  30.47200E-03 


(A) 


L-2000 


L-3000 


Retransmission 


Initial  Reflection 


-5.0000E+05 
-1  0000E+06 
-1  5000E+06 
-2.0000E+06 
-2.5000E+06 
-3.0000E+06 
-4  OOOOE+06 
-S.0000E+06 
-1  OOOOE+07 
-1  5000E+07 
-2  OOOOE+07 

Min :  -2.1 369E+07 
at  Node  19793 
Max:  1 .0601 E+05 
at  Node  11280 

P  (pa) 


Min :  5.0000E+03 
at  Cell  126,  1, 
Max:  1  4064E+06 
at  Cell  302,  1 .  2 


Figure  14. 


Pressure  Distribution  at 
Perfectly  Reflected  Boundary 
Solid  Bottom  Boundary  with 


t=30  msec  for:  (A) 
&  (B)  Lagrangian 
Retransmission 
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(A)  t=16  msec 


(B)  t=21  msec 


Stress-Z 


-5.0000E+05 
-1  0000E+06 
.1  5000E+06 
-2  OOOOE+06 
-2  5000E+06 
-3  0000E+06 
-4  OOOOE+06 
-5.0000E+06 
-1  0000E+07 
-1  5000E+07 
-2.0000E+07 

Mrt  -2  2677E-*07 
at  Node  19788 
Max:  7.5748E+06 
at  Node  3619 


Figure  15. 


Timeline  of  Vertical  (Z)  Stress  through 
Center  of  Lagrangian  Solid  Bottom  with  Nodal 

Constraints 


As  previously  discussed,  this  wave  is  the 
retransmission  wave.  A  review  of  the  vertical  (Z)  stress 
field  in  the  Lagrangian  solid  bottom  in  Figure  15  gives  a 
better  understanding  of  this  effect.  The  impact  of  the 
initial  shockwave  (Figure  15 (A) )  causes  an  increase  in 
compressive  stress  (blue  region) .  As  the  wave  travels 
through  the  material  it  is  reflected  off  the  nodal 
constraints  (Figure  15  (B) )  back  up  thru  the  model  causing  a 
decrease  in  the  compressive  stress  (pink  region) .  The 
expanded  portion  of  the  structure  then  causes  the 
retransmission  of  the  initial  shockwave  back  into  the  fluid 
domain.  The  Lagrangian  solid  bottom  then  continues  to 
compress  and  expand  as  it  attempts  dissipate  the  shock 
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energy.  In  nature,  once  the  shockwave  had  entered  the  soil 
it  would  propagate  through  the  infinite  soil  medium  until 
the  wave  dissipated  or  reached  a  different  medium,  like 
bedrock.  Due  to  the  finite  dimensions  of  the  model,  as  the 
shockwave  reaches  the  boundaries,  it  is  reflected  off  of 
the  fixed  nodal  constraints  and  causes  a  retransmission  to 
the  fluid  domain.  Simply  removing  the  nodal  constraints  on 
the  bottom  surface  will  only  serve  to  make  the  bottom  side 
of  the  structure  behave  as  a  free  surface  and  retransmit  a 
rarefaction  wave. 

2 .  Reducing  Retransmission 

Figure  14  (B)  makes  it  apparent  that  the  pressure  of 
the  retransmission  is  less  than  that  of  the  initial 

reflection.  This  is  a  result  of  energy  dissipation  as  the 
wave  passes  through  the  thickness  of  the  model.  Therefore 

to  minimize  the  retransmission  the  initial  solution  was  to 
increase  the  distance  through  which  the  wave  must  travel. 
By  increasing  the  bottom  thickness  the  retransmission  could 
be  reduced  to  a  negligible  level. 

A  follow-on  simulation  was  run  in  which  the  bottom 
thickness  was  increased  from  five  to  ten  meters.  The 
remainder  of  the  simulation  input  parameters  were  held 
constant  with  Table  10.  The  pressure-time  histories  for 
both  simulations  at  a  point  five  meters  above  the  bottom 

are  in  Figure  16.  By  doubling  the  thickness,  the 

retransmission  was  delayed  by  three  milliseconds  and 
decreased  in  pressure  by  2.5  x  105  Pa.  In  order  for  the 
retransmission  to  be  considered  negligible,  the  wave 
pressure  should  be  nearly  equal  to  the  hydrostatic  pressure 
the  bottom  would  normally  feel,  which  is  approximately  4.5 
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x  105  Pa.  Assuming  a  linear  relationship  between  bottom 
thickness  and  pressure  decrease,  the  bottom  would  have  to 
be  20  meters  thick  for  the  retransmission  pressure  to  be 
considered  negligible.  This  relationship  only  holds  for  a 
601b  charge  of  HBX-1.  An  increase  in  charge  size  would 
result  in  a  corresponding  increase  in  bottom  thickness.  The 
extraordinary  number  of  additional  structural  elements 
required  to  increase  the  thickness  would  require  a 
significant  increase  in  computational  resources  and  time. 


Figure  16.  Pressure-Time  History  5  Meters  above  Bottom 

(x=-5cm,  y=-5cm,  z=-3022cm) 

3 .  Non-Reflecting  Boundary  Segments 

The  boundary  condition  solution  had  to  not  only  keep 

the  Lagrangian  solid  bottom  in  position,  but  also  allow  the 

structure  to  behave  as  a  semi-infinite  domain.  The  problem 

of  modeling  a  semi-infinite  soil  medium  has  been 

encountered  previously  in  structural  dynamic  analyses  on 

dams  and  free-standing  structures.  The  approach  used  by 

O'Shea  utilized  LS-DYNA' s  transmitting  boundary  segments  in 
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order  to  simulate  a  semi-infinite  soil  domain  around  the 
Parker  Dam.  Dyna_N  contains  an  identical  function  listed  as 
non-reflecting  boundary  (NRB)  segments  [26]. 

According  to  the  Dyna_N  User's  Guide,  the  NRB  segments 
are  to  be  applied  to  the  exterior  boundaries  of  infinite 
domains  in  order  to  prevent  artificial  stress  wave 
reflections  generated  at  the  model  boundaries  from  re¬ 
entering  the  model  and  contaminating  the  results.  This  is 
done  through  the  use  of  impedance  matching  functions  that 
apply  normal  and  shear  stresses  of  the  form: 

^ normal  P^d^normal  &  ® shear  P^s^ tan 

r 

which  utilize  the  material  density  (p)  ,  the  dilatational 
wave  speed  (Cd)  ,  and  the  shear  wave  speed  (cs)  .  This  makes 
the  magnitude  of  these  stresses  proportional  to  the  normal 
and  tangential  particle  velocities  at  the  boundary  [27] . 

The  application  of  NRB  segments  allows  the  passage  of 
shockwaves  out  of  the  material.  Unfortunately,  due  to  their 
nature  as  Neumman,  or  flux,  boundary  conditions  they 
directly  conflict  with  the  application  of  fixed  nodal 
constraints.  Nodal  displacements  constraints  are  set 
through  enforcement  of  zero  acceleration  of  the  nodes  [19] . 
Since  the  applied  nodal  forces  are  simply  mass  multiplied 
by  acceleration,  the  application  of  both  boundary 
conditions  simultaneously  is  not  feasible.  It  is  unclear 
from  the  theoretical  manuals  for  Dyna_N,  which  boundary 
condition  is  given  precedent  when  both  conditions  are 
applied.  The  research  by  Zhenxia  proved  that  the 
transmitting  boundary  segments  in  LS-DYNA  are  insufficient 
to  handle  low-frequency  loading,  including  static  loads. 
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The  result  of  applying  solely  NRB  segments  to  a  static 
problem  is  a  permanent  translation  of  the  structure.  This 
result  is  confirmed  by  analysis  of  the  NRB  segments 


governing 

equations . 

When 

the 

object  is 

at 

rest. 

the 

particle 

velocities 

are 

zero. 

thus  the 

flux 

boundary 

conditions 

exert  zero 

force  on 

the  object. 

Not 

until 

the 

object  has  displaced  will  the  NRB  segments  fully  support 
the  static  loading  [28]  . 

4 .  Application  of  Nodal  Forces 

O'Shea  overcame  this  error  through  the  application  of 
static  forces  to  the  structure.  Initially,  the  model  was 
statically  loaded  with  the  semi-infinite  edges  set  as  fully 
constrained  nodes.  A  static  solution  was  found  for  the 
reaction  force  at  each  constrained  node.  These  static 
forces  were  then  applied  to  the  dynamic  model  to  act  as  the 
static  nodal  constraints.  Thus,  O'Shea  applied  two  flux 
boundary  conditions  to  his  structure:  the  NRB  segments  for 
the  shockwaves,  and  static  nodal  forces  to  replace  the 
nodal  constraints  [26]. 


Hydrostatic  Pressure 


\f 

\ . / 

Infinite  Bottom 


Soil  Structure  Free  Body  Diagram 
44 


Figure  17 . 


Table  11.  Soil  Nodal  Forces 


Nodal  Forces 
(x  1010  dynes) 

Interior 

Node 

Edge 

Node 

Corner 

Node 

F_z 

5.4130 

2.7065 

5.4130 

F_xy_l 

2.2489 

1.1244 

n/a 

F_xy_2 

4.6285 

2.3142 

n/a 

F_xy_3 

4.8246 

2.4123 

n/a 

F_xy_4 

5.0207 

2.5104 

n/a 

F_xy_5 

5.2169 

2.6084 

n/a 

F_xy_6 

2.6738 

1.3369 

n/a 

Figure  18.  Diagram  of  Lateral  Node  Force  Locations 

associated  with  Table  11 


A  simulation  was  conducted  using  a  five  meter  thick 
Lagrangian  solid  bottom  with  the  application  of  the  nodal 
force  method  used  by  O'Shea.  The  parameters  of  the 
simulation  are  located  in  Table  12.  In  addition  to  NRB 
segments  on  the  lateral  and  bottom  sides,  the  nodal  forces 
from  Table  11  were  applied  to  the  model.  These  nodal  forces 
were  calculated  based  upon  the  application  of  in-situ  soil 
hydrostatic  pressures  depicted  in  Figure  17.  Figure  19 
demonstrates  that  this  combination  of  boundary  conditions 
showed  almost  no  signs  of  retransmission  as  compared  with 
the  nodal  constraints  simulation.  Based  upon  these  results. 
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the  combination  of  NRB  segments  and  nodal  forces  will  be 
utilized  in  all  subsequent  models  investigating  additional 
bottom  modeling  parameters. 


Table  12.  Input  Parameters  for  Run  ID  4-07 


Charge 

Fluid 

Target 

Type 

HBX-1 

X  Width 

84  m 

Position 

N/A 

Weight 

60  Ibf 

Y  Width 

84  m 

Separation 

N/A 

Depth 

10  m 

Z  Depth 

35  m 

Bottom 

Dimensions 

Material 

Shape 

Parallelepiped 

Type 

Elastic 

X  Width 

84.02  m 

Density 

2.0  gram/cm3 

Y  Width 

84.02  m 

Elastic  Modulus 

2  x  10n  dyne/cm2 

Thickness 

5  m 

Poisson's  Ratio 

0.3 

#  Elements 

35,280 

EOS  Type 

None 

Boundary  Conditions 

All  nodes  on  lateral  faces  and  underside  of  the  Lagrangian  solid  bottom 
have  nodal  forces  normal  to  its  surface,  which  are  the  reaction  forces  of 
the  static  solution  of  the  structure.  The  same  faces  have  NRBs  applied. 

Figure  19.  Pressure  Time  History  5  meters  above  Ocean 

Bottom,  Euler  Coordinates  (x=-5cm,  y=-5cm,  z=- 
3022cm)  for  RunID_4_01  &  _07  Highlighting 
Effectiveness  of  NRB 
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The  initial  method  used  to  determine  the  nodal  forces 
was  a  rough  estimate  for  a  simple  geometric  shape  in  order 
to  validate  the  use  of  NRB  segments  and  nodal  forces.  In 
order  to  accurately  determine  the  true  nodal  forces,  a 
static  solution  of  the  structure  needed  to  be  solved.  This 
was  accomplished  through  the  use  of  the  ANSYS  12.1 
Structural  APDL  Software  Package.  For  the  nodal  forces  to 
be  easily  input  into  Dyna_N  it  was  critical  that  ANSYS  use 
the  exact  same  finite  element  model  that  was  generated  for 
Dyna_N.  This  required  converting  the  node  and  element  files 
generated  by  the  DYSMAS  processor  into  a  format  that  ANSYS 
could  read.  Conversely,  the  reaction  forces  for  the  nodal 
constraints  had  to  be  converted  into  the  required  Dyna_N 
format.  Both  of  these  tasks  were  accomplished  through  the 
development  and  application  of  MATLAB  script  files.  MATLAB 
was  utilized  due  to  the  author's  familiarity  with  it, 
though  this  perhaps  not  the  only  means  to  convert  and 
transfer  the  data  and  results.  Additionally  the  MATLAB 
script  generated  the  necessary  script  files  to  run  the 
ANSYS  solution  with  minimal  input  from  the  user.  Once  the 
node  and  element  file  had  been  read  into  ANSYS,  the 
hydrostatic  pressure  curve  was  generated  and  then  applied 
to  the  interface  surfaces.  The  script  then  constrained  the 
remaining  boundary  nodes  in  the  direction  normal  to  the 
surface  on  which  they  were  located.  For  example,  if  the 
node  lay  on  a  lateral  face  of  the  structure  that  was  normal 
to  the  X  axis,  then  ANSYS  constrained  that  node  from  moving 
in  the  X  direction.  ANSYS  then  determined  the  static 
solution  and  displayed  the  reaction  forces  at  the 
constrained  nodes. 


47 


A  note  of  caution  is  attached  to  the  use  of  NRB 
segments.  The  viscous  equations  the  NRB  segments  use  assume 
that  the  structure  is  composed  of  an  isotropic,  linear 
elastic  material.  All  simulations  to  this  point  have 
utilized  a  perfectly  elastic  material.  This  limitation 
could  prove  problematic  depending  upon  the  selection  of  the 
most  accurate  material  model  for  soil  [27] . 

C.  INITIAL  BOTTOM  WAVE  ELIMINATION 

Analysis  of  all  of  the  previous  simulations  shows  that 
at  the  start  of  the  simulation  a  pressure  wave  emanates  at 
the  water-soil  interface.  This  wave  formation  is 
demonstrated  in  Figure  20.  Figure  20(A)  is  the  state  of  the 
simulation  at  the  start  time.  Figure  20(B)  shows  the 
initial  wave  formation  with  the  pressure  dropping  0.5  x  10s 
Pa  at  the  wave  crest  within  the  first  1.5  milliseconds  of 
the  simulation.  After  another  three  milliseconds,  the 
pressure  has  increased  to  0.5  x  105  Pa  above  the  expected 
hydrostatic  value  in  Figure  20  (C)  .  The  variations  in  the 
contour  plots  in  Figures  20 (A-C)  demonstrate  that  this  wave 
is  present  along  the  entire  water-soil  interface.  Just 
before  the  initial  explosive  shockwave  interacts  with  this 
bottom  wave,  the  explosive  shockwave's  pressure  is 
approximately  3  x  106  Pa.  With  the  pressure  differential 
across  the  bottom  wave  of  a  third  the  magnitude  of  the 
explosive  shockwave,  this  flow  field  cannot  be  considered 
negligible  and  must  be  corrected. 
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(A)  t=  6.6  msec 


Time  :  I.U40M-M 


P  (pa) 


1  5000E+07 
4  6000E+05 
4.5000E+05 
4.4000E+05 
4.3000E+05 

4  2000E+05 
4.1000E+0S 
4.0000E+0S 
3  9000E+05 
3  8000E+05 
3.5000E+0S 
1  OOOOE+OS 

5  2000E+03 
0  0000E+00 


(B)  t=  8.0  msec 


(C)  t=  11.0  msec 
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1100 

■5 

&  1600 

o 

-2100 
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-3100 
-3600 

O.E+OO  l.E+05  2.E+05  3.E+05  4.E+05  5.E+05 

Pressure  (Pa) 


Min :  5  .0000E+03 
at  Cell  121,  1,2 
Max:  1 .0825E+07 
at  Cell  208,  1,47 


Figure  20. 


Initial  Bottom  Wave  Formation 


This 

structural 
simulation 
bottom  is 
pressure 
developed 


initial  bottom  wave  is  the  result  of  the 
deformation  of  the  interface.  When  the 
starts,  the  previously  unloaded  Lagrangian  solid 
instantaneously  loaded  with  the  hydrostatic 
from  the  fluid  domain  and  the  nodal  forces 
in  Section  IV. B. 4  to  act  as  nodal  constraints. 
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The  compressive  loading  of  the  structure  in  all  three 
dimensions  causes  the  structure  to  volumetrically  strain. 
This  strain  is  the  cause  of  the  displacement  of  the  water- 
soil  interface.  As  the  interface  displaces  downward,  the 
water  is  "pulled"  down  as  well  in  order  to  maintain  the 
interface.  Thus  the  initial  bottom  wave  takes  the  form  of  a 
tensile  pressure  wave.  As  with  any  elastic  material,  the 
structure  does  not  simply  deform  directly  to  its 
equilibrium  position.  Instead,  it  will  oscillate  like 
spring  about  its  equilibrium  position  until  it  settles  out 
[29]  .  This  oscillation  of  the  interface  in  the  vertical 
direction  causes  the  subsequent  pressure  rise  above 
hydrostatic  after  the  initial  tensile  wave.  This  motion  was 
observed  by  tracking  the  vertical  position  of  a  node  in  the 
Lagrangian  solid  bottom,  which  lay  on  the  water-soil 
interface  (Figure  21)  .  Points  A,  B,  and  C  in  Figure  21 
correspond  to  the  same  times  as  those  in  Figure  20.  An 
additional  factor  that  contributes  to  the  pressure  increase 
at  the  water-soil  interface  is  the  presence  of  back 
pressure  on  the  singly-wetted  interfaces  by  the  Gemini 
solver . 


Time  (sec) 


Displacement  of  Water-Soil  Interface 
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Figure  2 1 . 


In  order  to  eliminate  the  initial  bottom  wave,  the 
Lagrangian  solid  bottom  must  be  in  a  deformed  state  in 
equilibrium  with  the  hydrostatic  pressure  at  the  water-soil 
interface  and  the  applied  nodal  forces  developed  in  Section 
IV. B. 4  when  the  simulation  begins.  The  following  sections 
discuss  the  multiple  options  the  user  has  in  accomplishing 
this  task. 

Two  of  the  methods  to  reduce  the  initial  bottom  wave 
require  the  restart  of  a  DYSMAS  simulation.  Thus  an 
understanding  of  the  abilities  and  limitations  of 
restarting  simulations  is  required.  Coupled  DYSMAS 
simulations  can  be  restarted  with  the  addition  of  two  lines 
of  code  to  the  Dyna_N  input  file.  During  a  restart  several 
features  of  a  given  model  can  be  modified  from  the  original 
simulation.  Among  the  modifiable  features  are  changes  in 
termination  time,  deletion  of  materials,  deletion  of 
elements,  and  modifications  to  the  translational  or 
rotational  boundary  conditions  of  nodes.  The  restart 
simulations  utilized  restart  dump  files  created  at  specific 
time  intervals  throughout  normal  simulations  and  at  the 
conclusion  of  various  other  Dyna_N  processes  [27] . 

1 .  Long-Time  Static  Simulations 

Since  the  goal  is  to  establish  an  equilibrium  state  of 
the  model,  one  possible  solution  is  to  allow  the  structure 
to  interface  with  a  static  water  column  for  a  significantly 
long  duration  of  time.  At  the  conclusion  of  this  static 
simulation,  the  restart  file  produced  by  Dyna_N  is  then 
used  as  the  starting  structural  model  in  the  transient 
analysis  with  the  desired  UNDEX  event. 
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a.  Instantaneous  Loading 

The  first  simulation  was  conducted  with  a 
Lagrangian  solid  bottom  and  fluid  domain  with  the 
parameters  in  Table  13.  Nodal  displacement  constraints  were 
used  instead  of  nodal  forces  to  prevent  possibly 
significant  structural  movement  over  the  long  duration  of 
the  simulation.  The  nature  of  the  Dyna_N  restart  capability 
allows  the  user  to  reinsert  the  forces  during  the  UNDEX 
simulation.  The  differences  in  the  fluid  domain  were  the 
lack  of  an  explosive  charge,  the  application  of  wall 
boundary  conditions  on  the  lateral  sides,  and  setting  the 
back  pressure  for  the  singly-wetted  interfaces  to  zero.  The 
wall  boundary  conditions  were  applied  after  previous 
simulations  with  free  boundary  conditions  and  a  static 
water  column  had  shown  the  water  level  falling  several 
meters  over  the  course  of  the  simulations.  The  rational  for 
setting  the  back  pressure  to  zero  is  thoroughly  discussed 
in  Section  IV.C.2.d.  This  simulation  was  run  out  to  2.5 
seconds  with  the  expectation  that  the  structure  would 
achieve  equilibrium  within  that  time  span. 

The  results  of  this  instantaneously-loaded,  long¬ 
time  simulation  are  displayed  in  Figure  22.  On  the  left 
side  of  Figure  22,  it  is  important  to  note  that  the  maximum 
vertical  deformation  occurred  at  the  water-soil  interface 
with  a  value  of  -0.0072  centimeters.  This  closely  agrees 
with  an  ANSYS  Structural  APDL  static  solution  for  the  same 
structure  and  loading  conditions  where  the  deformation  of 
the  top  surface  was  -0.007  centimeters.  The  right  hand  side 
of  Figure  22  demonstrates  the  layered  nature  of  the 
vertical  stresses.  The  stresses  decrease  in  magnitude  going 
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from  the  interface  to  the  bottom  of  the  structure.  This 
behavior  is  consistent  with  the  uniaxial  compression  of  an 
object  that  is  prevented  from  expanding  in  the  other  two 
principal  directions.  The  ANSYS  solution  concurred  with 
this  stress  pattern  as  well.  While  this  method  appears  to 
be  reasonably  accurate,  a  quick  examination  of  the  vertical 
velocity  of  the  water-soil  interface  in  Figure  23  shows 
that  small  oscillations  are  still  occurring.  This  will 
cause  some  small  distortions  to  the  fluid  field. 

Table  13.  Run  ID  3-01  and  3-02  Simulation  Parameters 


Charge 

Fluid 

Target 

Type 

N/A 

X  Width 

84  m 

Position 

N/A 

Weight 

N/A 

Y  Width 

84  m 

Separation 

N/A 

Depth 

N/A 

Z  Depth 

35  m 

Bottom 

Dimensions 

Material 

Shape 

Parallelepiped 

Type 

Elastic 

X  Width 

84.02  m 

Density 

2.0  gram/cm3 

Y  Width 

84.02  m 

Elastic  Modulus 

2  x  10n  dyne/cm2 

Thickness 

5  m 

Poisson's  Ratio 

0.3 

#  Elements 

35,280 

EOS  Type 

None 

Boundary  Conditions 

All  nodes  on  lateral  faces  and  underside  of  the  Lagrangian  solid  bottom 
have  displacement  constraints  normal  to  the  surface  on  which  they  lie. 

NRB  segments  were  applied  to  the  same  faces. 
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Displacement-Z 


OOOOOE+OO 
-9.0231  E-04 
-1  8046E-03 
-2  7069E-03 
-3  6092E-03 
-4.5116E-03 
-5.41 39E-03 
-6.31 62E-03 
-7.2185E-03 


RunID  3-01:  Static  Run  to  2.5  sec 
Left:  Vertical  Deformation  (cm) 
Right:  Vertical  Stress  (dyne/cmA2) 


Stress-Z 


-3.3167E+06 
-3.4672E+06 
-3  6177E+06 
-3  7682E+06 
-3.91 87E+06 
-4.0692E+06 
-4  2197E+06 
-4.3702E+06 
-4  5208E+06 


Min  -7.21 85E-03 
at  Node  36298 


Min  -4.5208E+06 
at  Node  36142 


Figure  22.  Instantaneously  Loaded,  Long-Time  Static 

Solution 


Time  (sec) 


Figure  23. 


Interface  Vertical  Velocity  for  Run  ID  3-01 


b .  Ramped  Loading 

The  second  long-time  simulation  was  conducted 
with  conditions  applied  to  minimize  the  frequency  and 
magnitude  of  the  oscillations  the  model  experiences  prior 
to  equilibrium.  This  was  accomplished  by  gradually 


was 
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increasing  the  gravity  and  hydrostatic  pressure  loads  on 
the  structure  over  the  first  1.25  seconds  of  the 
simulation.  Implementing  this  effect  in  the  pre-processor 
was  relatively  straightforward.  The  load  curve  for  the 
gravitational  acceleration  was  modified  from  its  original 
constant  value  to  a  curve,  which  began  at  zero  magnitude  at 
time  zero  and  increased  linearly  to  its  full  value  at  1.25 
seconds  where  it  remained.  The  gradual  increase  in  pressure 
required  a  slightly  different  approach  because  the 
hydrostatic  loading  from  the  water-soil  interface  could  not 
be  directly  modified.  A  pressure  loading  in  addition  to  the 
hydrostatic  pressure  was  applied  to  the  interface  surface 
in  the  pre-processor.  This  pressure  was  set  to  be  of  equal 
magnitude,  but  opposite  direction  to  the  hydrostatic 
loading.  This  pressure  loading  was  given  a  load  curve  that 
started  at  full  magnitude  at  time  zero  and  decreased 
linearly  to  zero  pressure  at  1.25  seconds  where  it 
remained . 

Ramping  the  loading  of  the  structure  demonstrated 
a  significant  increase  in  accuracy.  The  maximum  deflection 
and  vertical  stress  layers  in  Figure  24  are  nearly 
identical  to  those  in  Figure  22.  The  true  benefit  of  the 
ramped  loading  is  the  final  velocity  of  the  interface 
surface  is  only  0.001  cm/ s  in  Figure  25  vice  the  -0.457 
cm/ s  at  the  end  of  the  instantaneously  loaded  structure  in 
Figure  23. 
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Figure  24.  Ramp  Loaded,  Long-Time  Static  Solution 


Time  (sec) 


Figure  25.  Interface  Vertical  Velocity  for  Run  ID  3-02 


c .  Summary 


Even  though  Dyna_N  is  not  the  most  suitable 
solution  for  static  problems,  the  method  of  long-time 
simulations  to  achieve  a  deformed  static  structure  is 
viable.  The  primary  advantage  of  this  method  is  that  there 
are  no  direct  conflicts  with  any  material  models  or 
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combination  of  structures.  The  obvious  drawback  is  that  is 
computationally  expensive  to  statically  simulate  large 
structural  models  over  a  long  time.  This  computational 
expense  led  to  the  examination  of  additional  technigues 
that  would  eliminate  the  initial  bottom  wave  [27] . 

2 .  Dynamic  Relaxation 

The  Dynamic  Relaxation  (DR)  tool  in  Dyna_N  and  the 
associated  Prestress  program  in  Gemini  were  specifically 
designed  for  the  problem  of  finding  the  deformed  shape  and 
the  associated  stresses  of  a  finite  element  model  in  static 
equilibrium  prior  to  the  start  of  an  UNDEX  simulation  [18], 
[27]  .  The  DR  algorithm  is  based  upon  the  fact  that  the 
solution  of  a  damped  dynamic  solution  converges  on  a  quasi¬ 
static  solution  as  time  approaches  infinity.  After  applying 
the  defined  static  and  interface  loads,  the  solver  steps 
through  the  algorithm  determining  the  kinetic  energy  of  the 
structure  and  tracking  the  maximum  total  kinetic  energy. 
Once  the  solver  determines  that  the  current  kinetic  energy 
of  the  system  is  a  user-defined  percentage  of  the  maximum 
total  kinetic  energy  the  system  has  seen,  the  solver  stops 
the  algorithm  and  writes  a  Dyna_N  restart  file  containing 
the  final  deformed  shape  and  internal  stresses  at  the  last 
solution  step.  The  default  percentage  of  maximum  kinetic 
energy  is  0.1  percent  [27]  .  In  order  for  a  Dyna_N  restart 
file  to  be  coupled  with  a  dynamic  fluid  simulation,  the 
restart  file  must  be  created  as  a  coupled  run.  For  DR 
solutions  this  coupling  is  accomplished  through  the 
application  of  the  Gemini  Prestress  function,  which 
provides  the  hydrostatic  pressure  loading  for  the  interface 
segments.  Unlike  a  normal  Gemini  fluid  simulation,  Gemini 
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Prestress  is  more  computationally  efficient  because  it  does 
not  create  and  solve  an  entire  fluid  domain  [18] . 


a.  Limitations 

The  use  of  DR  and  Gemini  Prestress  is  not  without 
limitations.  The  DR  algorithm  is  incompatible  with  non-zero 
displacement  or  any  velocity  boundary  conditions.  This  is 
in  direct  conflict  with  the  application  of  NRB  segments,  as 
they  are  velocity-based  boundary  conditions.  Additionally, 
the  nature  of  restart  files  does  not  allow  the  addition  of 
NRB  segments  after  the  DR  solution.  This  conflict  was 
proven  valid  after  no  DR  simulation  could  be  started  when 
the  structure  contained  NRB  segments.  The  solver  returned 
an  error  of  excessive  deformation  or  improper  definition  of 
a  solid  element.  As  the  material  used  was  a  perfectly 
elastic  material  it  is  not  possible  to  have  excessive 
deformation  [19]. 

A  second  limitation  noted  in  the  literature  is 
the  potential  for  overshoot.  If  the  loads  are  applied  too 
rapidly,  the  structure  may  oscillate  about  the  solution 
several  times  prior  to  achieving  the  minimum  kinetic 
energy.  This  can  have  a  severe  impact  upon  materials  that 
are  history-dependent.  Soils  that  can  experience  compaction 
are  prime  examples  of  history-dependent  materials.  While 
this  research  solely  utilizes  perfectly  elastic  material, 
this  fact  is  important  to  understand  prior  to  the 
application  of  more  accurate  soil  material  models.  The 
result  of  the  overshoot  would  be  that  the  soil  is  over¬ 
compacted  prior  to  the  UNDEX  simulation  [27] . 

Several  Lagrangian  solid  bottoms  with  varying 

boundary  conditions  were  dynamically  relaxed.  The  results, 
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in  some  cases  were  promising,  but  ultimately  proved 
unfeasible  for  application  to  Lagrangian  solid  bottoms.  The 
following  is  a  discussion  of  the  parameters  and  results  of 
each  variation  attempted. 

Jb.  With  Nodal  Forces 

A  DR  simulation  utilizing  only  static  nodal 
forces  at  the  boundaries  was  computationally  achieved,  but 
physically  inaccurate.  The  lack  of  nodal  displacement 
constraints  allowed  the  structure  to  translate  vertically 
nearly  60  centimeters  before  a  solution  was  converged  upon 
in  Figure  26.  Further  examination  of  the  results  showed 
that  the  model  was  not  at  rest,  but  moving  vertically  with 
a  velocity  of  four  centimeters  per  second. 


Displacement-2 


5.9984E+01 
5.9954E+01 
5.9924E+01 
5.9894E+01 
5.9864E+01 
5.9835E+01 
5.9805E+01 
5.9775E+01 
5.9745E+01 

Min :  5.9745E+01 
at  Node  39740 
Max:  5.9984E+01 
at  Node  1 


RunlO_3_03  DR 

BC:  Static  Node  Forces  Only 

Of  Note:  Minimum  Vertical  Displacement  60  cm 


Figure  26.  DR  with  Static  Nodal  Forces  Only 
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c.  With  Displacement  Constraints 

While  the  application  of  nodal  forces  failed,  the 
ability  to  modify  nodal  forces  and  constraints  in  a  restart 
file,  allows  the  user  to  dynamically  relax  a  structure  with 
displacement  constraints  and  then  replace  them  with  nodal 
forces  in  the  restart  of  the  transient  analysis.  The 
displacement  constraints  were  applied  on  the  lateral  and 
underside  of  the  Lagrangian  solid  bottom.  The  nodal 
constraints  prevented  movement  in  the  direction  normal  to 
the  surface  where  the  node  was  located.  Once  again,  a  DR 
solution  was  computationally  achieved,  but  physically 
inaccurate.  Observing  the  vertical  displacement  through  the 
structure  in  Figure  27,  the  bottom  surface  is  fixed  at  zero 
displacement  as  expected.  The  inaccuracy  is  noted  in  the 
near-zero  displacement  of  the  top  surface,  which  should 
have  the  greatest  displacement.  After  discussions  with  a 
DYSMAS  subject  matter  expert,  it  was  determined  that  this 
was  likely  due  to  the  back  pressure  applied  to  the  singly- 
wetted  interfaces  (SWI)  by  the  Prestress  software.  During  a 
transient  analysis,  the  user  has  the  ability  to  adjust  the 
back  pressure.  The  Gemini  Prestress  function  does  not  give 
the  user  this  option.  The  back  pressure  defaults  to  1  x  106 
dyne/ cm2  [ 30 ] . 
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RunID  3-03  DR 

XY  and  Z_min  Faces  Constrained  in  Normal  Directions 
Of  Note:  Interior  Nodes  Displace  Further  than  Top  Surface 


Min :  -1 .3250E-05 
at  Node  5 
Max:  8.6537E-04 
at  Node  28905 


Figure  27 


DR  with  Nodal  Constraints 


d.  Application  of  Counter  Pressure 

The  final  DR  simulation  utilized  the  nodal 
displacement  constraints  with  the  addition  of  a  pressure 
loading  to  counteract  the  effects  of  the  Gemini  Prestress 
back  pressure.  The  counter  pressure  was  applied  through  the 
application  of  compressive  pressure  loading  on  the  singly- 
wetted  interface  of  1  x  106  dyne/cm2.  A  profile  of  the 
vertical  deformation  through  the  thickness  of  the  model  is 
shown  in  Figure  28.  In  this  case,  the  interface  surface 
displaced  the  most  at  -0.0018  centimeters.  The  ANSYS  static 
solution  determined  that  the  displacement  of  the  surface 
should  be  -0.007  centimeters.  This  set  of  boundary 
conditions  provided  the  closest  representation  of  the 
static  structure. 
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Figure  28.  DR  with  Nodal  Constraints  and  Counter 

Pressure 


e .  Summary 

While  the  application  of  nodal  displacement 
constraints  and  counter  pressure  gives  a  reasonable  DR 

solution,  the  conflict  with  the  NRB  segments  significantly 
diminishes  the  benefits  of  using  dynamic  relaxation  in 

Lagrangian  solid  bottom  modeling.  An  even  greater  failure 
of  dynamic  relaxation  is  its  apparent  inability  to 

successfully  relax  two  separate  structures,  such  as  the 
ocean  floor  and  the  FSP.  A  DR  simulation  containing  two 

structures  failed  to  converge  on  a  solution.  While  the  runs 
experienced  no  errors,  the  kinetic  energy  requirement  was 
simply  never  met.  Dynamic  relaxation  appears  to  be  a  semi- 
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effective  tool  when  modeling  the  ship  structure  or  the 
bottom,  not  both  simultaneously. 

3.  Dyna_N  Prestress 

The  application  of  the  Dyna_N  Prestress  function 

offered  the  best  method  to  reduce  the  effect  of  the  initial 
bottom  wave.  The  Prestress/Prestrain  function  in  Dyna_N 

allows  the  user  to  define  the  stress  tensor  (oxx,  oyy,  ozz, 
oxy ,  ayz,  ozx)  and  the  equivalent  plastic  strain  (psv)  at  the 
simulation  start  time  for  any  element  in  the  structure 
[27] .  This  function  is  only  briefly  mentioned  in  the  Dyna_N 
User' s  Manual  and  the  pre-processor  does  not  have  an  input 
interface  for  the  user  to  utilize  the  function.  After 
discussions  with  subject  matter  experts,  this  function  was 
added  to  the  code  by  the  DYSMAS/P  co-developer  IABG,  and 

not  the  Dyna_N  code  owners  at  Lawrence  Livermore  National 

Laboratory.  IABG  designed  the  Dyna_N  Prestress  function  to 
assist  in  simulations  involving  metal  forming.  As  it  was 

designed  for  use  with  metal,  the  function  only  works 

correctly  when  applied  to  elastic  or  elastic-plastic 
material  models  [31]  .  While  the  exact  algorithm  was  not 
available  for  reference,  it  is  theorized  that  the  stress 
tensor  serves  as  the  starting  point  for  the  elemental 

stress  calculations  in  Dyna_N  at  the  first  step.  It  does 
not  appear  that  this  stress  tensor  is  translated  into  an 
initial  elastic  strain  tensor. 

The  initial  stress  tensors  utilized  were  initially 
calculated  by  hand  in  order  to  test  the  application  of  the 
Dyna_N  Prestress  function.  The  stress  tensors  were  made  on 
the  assumption  that  the  only  stress  in  the  soil  was 

hydrostatic.  With  the  generic  flat  Lagrangian  solid  bottom, 
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which  was  five  meters  thick  with  one  meter  thick  elements, 
a  simple  pressure  calculation  was  done  at  the  center  of 
each  element  to  determine  the  pressures  of  each  row  in 
Table  14.  The  input  card  for  the  Dyna_N  Prestress  function 
was  generated  using  a  MATLAB  script.  All  shear  portions  of 
the  stress  tensor  were  assumed  to  be  zero. 


Table  14.  Bottom  Hydrostatic  Pressure  Distribution 


Layer 

Depth  From 
Air-Water 

Interface 

(cm) 

Depth  From 
Water-Soil 

Interface 

(cm) 

Pressure  At 
Min.  Depth 
(dyne/cm2) 

Pressure  At 
Max.  Depth 
(dyne/cm2) 

Average  Layer 
Pressure 
(dyne/cm2) 

1 

3500-3600 

0-100 

4.4323  x  106 

4.6285  x  106 

4.5303  xlO6 

2 

3600-3700 

100-200 

4.6285  x  106 

4.8246  x  106 

4.7265  xlO6 

3 

3700-3800 

200-300 

4.8246  x  106 

5.0207  x  106 

4.9227  x  106 

4 

3800-3900 

300-400 

5.0207  x  106 

5.2169  x  106 

5.1188 x 106 

5 

3900-4000 

400-500 

5.2169  x  106 

5.4130  x  106 

5.3149  x  106 

The  average  layer  pressure  in  Table  14  was  used  as  the 
normal  stress  for  the  Dyna_N  Prestress  inputs.  A  prestresed 
simulation  was  conducted  with  all  of  the  same  parameters  as 
the  simulation  in  Table  13.  The  results  in  Figure  29  were 
promising.  The  application  of  prestress  significantly 
reduced  the  initial  pressure  drop.  Even  though  the 
prestress  function  does  not  deform  the  structure  prior  to 
problem  start,  it  does  appear  to  stiffen  the  structure. 
This  added  stiffness  allows  the  structure  to  settle  at  a 
slower  rate,  which  diminishes  the  initial  pressure  drop. 

Even  with  the  prestress,  the  bottom  still  emits  a  wave 
at  simulation  start.  This  is  due  to  the  use  of  simple 
hydrostatic  calculations  to  determine  the  prestress  tensor 
and  the  presence  of  back  pressure  on  the  singly-wetted 
interface  (SWI) .  In  order  to  more  accurately  represent  the 
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prestress  tensor,  the  static  solution  in  ANSYS  was  once 
again  utilized.  The  stress  tensor  for  every  element  in  the 
structure  was  printed  from  the  ANSYS  static  solution.  A 
MATLAB  script  was  used  to  convert  the  ANSYS  output  into  the 
format  of  the  Dyna_N  Prestress  input.  Follow-on  simulations 
that  utilized  the  ANSYS-generated  Dyna_N  Prestress  inputs 
yielded  results  consistent  with  Figure  29(B). 


(A)  —  No  Prestress 

(B)  -  Prestressed 
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Figure  29.  Effect  of  Dyna_N  Prestress  on  Bottom  Wave 

4 .  Back  Pressure  Elimination 

The  remaining  portion  of  the  initial  bottom  wave  in 
Figure  29(B)  is  due  to  the  presence  of  back  pressure  on  the 
singly-wetted  interface.  The  back  pressure  is  skewing  the 
transient  results  in  the  same  fashion  in  which  it  distorted 
the  dynamic  relaxation  solution.  There  are  two  solutions  to 
this  problem.  The  simplest  to  apply  is  the  elimination  of 
back  pressure  in  the  Gemini  input.  This  only  requires  the 
modification  of  a  single  line  of  code.  A  drawback  occurs 
with  the  addition  of  a  ship  model.  The  FSP  model  was  turned 
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into  a  singly-wetted  interface  by  adding  an  artificial 
roof.  This  freed  the  user  from  establishing  interfaces  in 
the  interior  of  the  model.  The  elimination  of  the  back 
pressure  would  necessitate  a  change  in  the  ship  model. 

A  second  method  to  eliminate  the  back  pressure  is  to 
apply  a  pressure  that  is  of  equal  magnitude  but  opposite 
direction  to  the  back  pressure.  This  is  done  through  the 
application  of  pressure  loading  on  the  interfaces,  which  is 
easily  accomplished  in  the  preprocessor.  Conveniently,  the 
pressure  loading  is  applied  on  a  separate  input  file  from 
the  nodal  forces.  This  ensures  there  is  no  input  conflict 
between  the  pressures  and  the  nodal  reaction  forces.  This 
method  of  back  pressure  elimination  was  applied  to  all 
further  simulations. 

5.  Benefits  of  Decreased  Bottom  Thickness 

The  final  recommendation  on  minimizing  the  bottom  wave 
is  to  decrease  the  thickness  of  the  model.  The  cause  of  the 
initial  bottom  wave  was  the  movement  of  the  top  surface  of 
the  Lagrangian  solid  bottom.  This  is  the  byproduct  of  the 
structural  volumetric  strain  due  to  the  initial  hydrostatic 
loading.  While  the  hydrostatic  pressure  will  always  control 
the  amount  of  strain,  the  range  of  the  displacement  of  the 
structure  can  be  minimized  by  decreasing  the  thickness  of 
the  structure.  Thus  the  magnitude  of  the  initial  bottom 
wave  is  diminished.  This  concept  was  validated  by  running 
two  sets  of  simulations  in  which  the  only  variable  was  the 
thickness  of  the  model.  The  first  set  used  all  of  the 
parameters  developed  in  this  chapter.  The  only  difference 
in  the  second  set  was  that  the  Dyna_N  Prestress  inputs  were 
not  included.  Figure  30  is  a  graph  of  the  pressure 
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deviation  from  hydrostatic  as  a  function  of  distance  from 
the  bottom  for  the  first  set  of  simulations.  In  both  cases, 
three  bottom  thickness  of  two,  five,  and  ten  meters  were 
applied.  Figure  31  graphs  the  same  results  for  the  set 
without  Dyna_N  Prestress  inputs.  In  all  cases  the 
individual  element  sizes  were  kept  constant  at  one  cubic 
meter . 


%  AP  (t=2m)  %  AP  (t=5m) - %AP(t=10m) 


Depth  (cm) 


Figure  30.  Effect  of  Bottom  Thickness  on  Initial  Bottom 

Wave  with  Dyna_N  Prestress 
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Figure  31 .  Effect  of  Bottom  Thickness  on  Initial  Bottom 

Wave  without  Dyna_N  Prestress 

From  Figure  30,  it  is  apparent  that  when  the  Dyna_N 
Prestress  inputs  are  provided  the  initial  bottom  wave  is 
only  0.1%  greater  than  hydrostatic  pressure.  Therefore,  the 
Lagrangian  solid  bottom  can  be  made  a  thin  as  possible  in 
order  to  minimize  the  number  of  elements  and  decrease  the 
simulation  time.  Figure  31  highlights  the  benefit  of 
minimizing  the  thickness  when  Dyna_N  Prestress  inputs 
cannot  be  provided.  As  the  thickness  decreases,  the 
relative  magnitude  of  the  initial  bottom  wave  is  decreased. 
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D .  RECOMMENDED  METHOD  FOR  LAGRANGIAN  SOLID  BOTTOM 

MODELING 

Constructing  the  Lagrangian  solid  bottom  outside  of 
the  fluid  domain,  as  diagrammed  in  Figure  13,  encompasses 
most  of  the  potential  hurdles  to  implementing  a  Lagrangian 
solid  bottom.  With  the  correct  application  of  initial  and 
boundary  conditions,  this  approach  offers  a  more 
computationally  efficient  method  of  Lagrangian  solid  bottom 
modeling  than  completely  containing  the  bottom  within  the 
fluid  domain. 

This  chapter  has  explored  several  variations  of 
initial  and  boundary  conditions  for  the  Lagrangian  solid 
bottom  model.  The  best  combination  of  these  parameters  is 
listed  in  Table  15.  The  significant  exception  to  this  study 
is  the  material  model  used  for  to  simulate  the  soil.  A 
perfectly  elastic  material  with  parameters  roughly  similar 
to  soil  was  utilized  in  order  to  minimize  material  model 
specific  errors.  Further  research  will  be  required  to 
determine  the  best  material  model  to  use.  The  remainder  of 
this  parametric  study  in  the  effects  of  an  explicitly 
modeled  structural  ocean  floor  on  a  shallow  water  UNDEX 
event  will  apply  these  parameters. 


69 


Table  15.  Best  Lagrangian  Solid  Bottom  Modeling  Parameters 


Parameter 

Best  Choice 

Lateral 

Dimension 

Just  greater  than  lateral  dimensions  of  fluid 
domain . 

Vertical 

Dimension 

Thinner  is  more  computationally  efficient  and 
required  if  use  of  Dyna  N  Prestress  inputs  is 
not  feasible. 

Fluid-Soil 

Vertical 

Overlap 

Overlap  must  be  large  enough  to  ensure  that 
the  vertical  deformation  of  the  Lagrangian 
solid  bottom  does  not  shift  the  water-soil 
interface  out  of  the  fluid  domain. 

Boundary 

Conditions 

Nodal  Forces  that  correspond  to  the  reaction 
forces  required  for  nodal  displacement 
constraints  on  surfaces  not  interfaced  with 
the  fluid. 

Non-reflecting  Boundary  Segments  on  surfaces 
not  interfaced  with  the  fluid  domain. 

Counter  Pressure  on  the  singly-wetted 
interfaces  is  necessary  if  the  ship  model 
utilizes  singly-wetted  interfaces  as  well. 

Initial 

Condition 

Dyna  N  Prestress  input  for  every  solid  element 
derived  from  the  element  stress  tensors  of  the 
static  solution,  if  and  only  if  the  material 
model  used  in  the  Lagrangian  solid  bottom  is 
elastic  or  elastic  plastic. 

Use  restart  file  from  a  long-time,  ramp  loaded 
static  simulation  if  the  material  model  is  not 
elastic  or  elastic-plastic. 
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V.  COMPARISON  OF  LAGRANGIAN  SOLID  BOTTOM  TO 
CURRENT  BOTTOM  MODELS 


A.  CURRENT  MODELING  METHODS  AND  LIMITATIONS 

There  are  two  currently  accepted  ways  in  which  to 
model  the  ocean  bottom  in  the  DYSMAS  software  suite.  All  of 
these  methods  are  implemented  through  the  Eulerian  fluid 
solver,  Gemini.  The  first  is  the  imposition  of  a  perfectly 
reflective  boundary  condition  on  the  fluid  domain.  The 
second  is  the  use  of  a  soil-like  material  layer  at  the 
bottom  of  the  fluid  domain.  Both  have  their  own  unique 

limitations . 

1 .  Wall  Boundary  Condition 

The  simplest  method  is  to  enforce  a  perfectly 
reflective,  wall  boundary  condition  on  the  bottom  face  of 
the  fluid  domain.  The  pre-process  for  this  method  is  a 
simple  one  line  input  in  the  Pregemini  input.  When 

conducting  simulations  involving  gravity,  a  wall  boundary 
condition  is  already  required  in  order  to  keep  the  water 
from  draining  out  of  the  bottom  of  the  fluid  domain.  The 
only  boundary  condition  adjustment  available  is  the 
modification  of  the  amount  of  the  reflection  from  zero  to 
100  percent  of  the  incident  pressure.  While  this  option  is 
the  easiest  to  implement,  the  limitations  are  obvious.  The 
wall  condition  acts  as  a  perfectly  rigid  boundary,  unlike  a 
soil  whose  response  to  pressure  loading  can  range  from 

elastic  to  plastic.  The  creation  of  a  three  dimensional 

Cartesian  fluid  domain  requires  that  the  bottom  surface  be 
flat.  This  prevents  this  boundary  condition  from  being  able 
to  form  contoured  bottom  structures  [18]. 
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2. 


Euler  EOS  Material  Layer 


The  bottom  modeling  method  which  is  currently  the 
accepted  standard  when  using  DYSMAS  in  a  shallow  water 
environment  utilizes  various  Euler  EOSs.  During  the 
creation  of  the  fluid  domain,  in  addition  to  air  and  water, 
a  layer  of  soil  is  created  at  the  bottom  of  the  fluid 
domain.  The  EOSs  for  this  soil  are  either  a  Mie-Griineisen 
or  a  P-alpha  form.  The  Mie-Griineisen  is  an  EOS  relating  the 
pressure,  energy,  and  density  of  the  material  through  the 
shock  speed  versus  particle  speed  profile  of  the  material. 
The  P-alpha  EOS  combines  the  Mie-Griineisen  EOS  with  an  air 
EOS  in  order  to  more  accurately  simulate  a  material  with 
collapsible  porosity  [18]  .  This  method  of  bottom  modeling 
has  been  validated  and  used  extensively  in  previous  UNDEX 
simulations  conducted  in  DYSMAS  to  date  [8],  [9],  [10], 
[11]  .  While  it  is  the  standard,  there  are  obvious 
drawbacks.  As  an  Eulerian  fluid,  the  soil  material  cannot 
support  shearing  forces,  whereas  actual  soil  does  support 
shear  loading.  A  new  version  of  the  Gemini  code  which 
allows  for  the  application  of  a  viscous  fluid  is  in 
development,  but  was  not  available  at  the  time  of  this 
research.  Creating  a  contoured  bottom  using  an  Euler  EOS  is 
not  accurate.  The  reason  for  this  limitation  is  explained 
further  in  Chapter  VI . 

B.  COMPARISON  OF  BOTTOM  MODELING  METHODS 

A  set  of  simulations  were  conducted  to  compare  the 
effects  that  each  accepted  bottom  modeling  method  and  the 
Lagrangian  solid  bottom  method  developed  in  Chapter  IV  had 
on  the  fluid  domain.  Each  simulation  used  the  standard 
charge  size  and  FSP  placement  location  from  Chapter  III. 
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The  first  simulation  utilized  a  Gemini  wall  boundary 
condition  with  a  100  percent  reflection.  The  second  bottom 
model  was  a  Mie-Griineisen  EOS  for  a  clay  soil.  This  soil 
layer  was  five  meters  thick,  with  a  wall  boundary  condition 
applied  on  the  underside  of  the  soil.  The  final  model  was 
generated  from  the  best  practices  developed  in  Chapter  IV. 
The  model  was  five  meters  thick  and  consisted  of  an  elastic 
material  model  whose  properties  were  consistent  with  a  clay 
soil  provided  by  Bangash  [23]  .  While  the  FSP  was  included 
in  each  simulation,  this  initial  comparison  was  solely 
focused  on  the  fluid  domain  response  to  each  model.  Each 
simulation  was  run  out  to  one  second. 

1.  Effect  on  Initial  Bottom  Reflection 

The  first  comparison  point  was  taken  at  35 
milliseconds  after  the  charge  was  detonated.  Figures  32 
thru  34  represent  the  pressure  distribution  through  the 
water  column  for  the  wall  condition,  Euler  soil,  and 
Lagrangian  solid  bottom  models  respectively.  This 
particular  time  was  chosen  because  it  captures  the  entire 
bottom  reflection  response  of  each  model  prior  to  the 
reflected  shockwave  impacting  the  bulk  cavitation  zone 
(seen  in  white)  .  The  black  line  across  the  bottom  of  the 
plot  in  Figure  33  is  interface  between  the  water  and  the 
Euler  soil.  A  similar  line  in  Figure  34  is  the  interface 
between  the  water  and  Lagrangian  solid  bottom.  Of 
significance  in  Figure  33  is  the  presence  of  a  double 
reflection.  The  pressure  wave  closest  to  the  surface  is  the 
true  bottom  reflection.  The  second  reflection  is  a  result 
of  the  wall  condition  placed  beneath  the  soil  layer.  In  all 
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cases 


the  timing  of  the  initial  bottom  reflection  appears 
to  be  consistent.  None  of  the  models  appears  to  delay  the 
reflection . 


RunID  1-00:  Bottom  Wall  Boundary  Condition 
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Figure  32.  Pressure  at  35  msec  for  a  wall  boundary 
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RunID  1-01:  Euler  Clay  Bottom 
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Figure  33. 
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Figure  34.  Pressure  at  35  msec  for  Solid  Bottom 
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Figure  35.  Bottom  Reflection  Pressure  Comparison 


A  more  accurate  comparison  was  drawn  by  plotting  the 
pressure  versus  depth  at  35  milliseconds  for  all  three 
cases  in  Figure  35.  Of  note  is  the  clear  correlation  in 
shape  between  the  wall  boundary  condition  and  the 
Lagrangian  solid  bottom.  The  only  difference  between  the 
two  appears  to  be  the  magnitude  of  the  reflected  pressures. 
The  Euler  soil  looks  to  follow  the  same  pattern  from  the 

surface  to  -1000  centimeters.  At  that  point  the  response 
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becomes  skewed  by  the  second  reflection.  The  influence  of 
the  second  reflection  can  be  reduced  to  negligible  values 
with  an  increase  in  thickness  of  the  Euler  clay  layer  [32] . 
This  reduction  would  come  at  the  cost  of  increased 
computational  resources.  Comparing  the  magnitudes  of  the 
initial  reflections,  the  wall  condition  provides  the 
greatest  reflection  and  the  Euler  soil  provides  the  least. 
The  elastic  material  model  used  in  the  Lagrangian  solid 
bottom  could  be  adjusted  in  order  to  reflect  a  pressure 
wave  of  equal  magnitude  to  that  of  the  Euler  soil. 
Replacing  the  elastic  material  with  a  more  representative 
soil  model  could  also  improve  upon  this  response. 

2 .  Effect  on  Gas  Bubble 

The  bottom  model  type  also  affected  the  bubble 
response.  Table  16  compares  the  response  of  each 
simulation'’ s  bubble  and  associated  first  pulse.  Once  again 
the  Lagrangian  solid  bottom  closely  follows  the  wall 
boundary  condition  in  all  areas.  The  significant  outlier  is 
the  incident  pressure  of  the  Euler  soil's  first  pulse, 
which  is  nearly  double  the  incident  pressure  of  the  other 
two  models. 


Table  16.  Max  Radius  and  First  Pulse  Comparison 


Model 

Max  Bubble 
Radius  (cm) 

1st  Pulse  on 
FSP  (msec) 

1st  Pulse 
Pressure  (Pa) 

Empirical 

418 

552 

N/A 

Wall  B.C. 

(Run  1-00) 

475 

563 

6  x  105 

Euler  Soil 
(Run  1-01) 

482 

576 

12  x  105 

Lag.  Solid  Bottom 
(Run  7-06) 

475 

563 

7  x  105 
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VI.  EFFECT  OF  BOTTOM  CONTOURS  ON  UNDEX  EVENT 


Modeling  the  bottom  as  a  structure  in  Dyna_N  provides 
the  capability  to  create  bottoms  which  are  not  flat.  When 
modeling  the  bottom  in  Gemini,  it  is  not  possible  to 
establish  three  dimensional  bottom  shapes  that  are  in 
equilibrium.  It  is  possible  to  fill  an  arbitrary  shape  in 
the  Euler  grid  with  a  soil  material.  Two  limitations 
prohibit  the  soil  from  staying  in  the  arbitrary  shape.  The 
first  is  the  density  mismatch  between  the  soil  and  water 
EOS.  Since  the  soil  is  treated  as  an  inviscid  Eulerian 
fluid,  over  time  the  denser  soil  will  settle  to  the  bottom 
of  the  fluid  geometry.  Several  simulations  without 
explosive  charges  were  run  to  determine  if  this  effect  was 
visible.  After  only  100  milliseconds,  the  soil  was  observed 
to  be  settling  in  Figure  36.  The  left  side  of  Figure  36  is 


the  density 

profile 

for  the 

domain.  The  right 

side 

of 

Figure  36 

is  the 

vertical 

velocity  profile 

at 

100 

milliseconds 

The  blue  region 

indicates  that 

the 

soil 

columns  are  moving  down  at  a  rate  of  nearly  two  centimeters 
per  second.  The  pink  region  is  water  moving  upward  at  two 
centimeters  per  second.  The  actual  displacement  over  the 
course  of  one  second  is  negligible  to  the  overall  geometry. 
The  root  cause  of  this  movement  is  an  uneven  pressure 
distribution  through  the  fluid  domain.  This  pressure 
distortion  has  severe  consequences  for  the  simulation  and 
cannot  be  avoided  when  using  Euler  soil. 
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Figure  36. 


Density  and  Vertical  Velocity  Profile  of 
Euler  Contoured  Bottom 


In  order  to  create  a  three  dimensional  bottom  in 
Gemini,  the  user  must  utilize  the  BODYFILL  command. 
BODYFILL  allows  the  user  to  fill  an  arbitrary  shape  in  the 
fluid  domain  with  a  given  material  [18] .  Once  the  shape  has 
been  filled,  Gemini  applies  gravitational  forces  to  the 
filled  material.  The  complication  is  that  Gemini  cannot 
equalize  the  pressure  in  the  fluid  domain  with  the  filled 
shape.  This  results  in  the  fluid  domain  pressure 
distribution  depicted  in  Figure  37.  The  pressure  gradient 
in  the  soil  is  much  greater  than  the  water  due  to  the 
difference  in  density.  The  pressure  mismatch  on  the 
vertical  faces  and  at  the  bottom  of  the  channel  creates 
considerably  large  pressure  waves  when  the  simulation 
begins.  Just  as  the  initial  bottom  wave  in  Section  IV.  C 
distorts  the  explosive  shockwave,  so  too  does  this  pressure 
mismatch.  Unlike  the  Lagrangian  solid  bottom,  there  is  no 
method  to  correct  this  distortion  in  the  Eulerian  solver. 
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Figure  37.  Pressure  Gradient  in  Euler  Contoured  Bottom 


In  Chapter  IV,  this  thesis  explored  the  various 
methods,  distortions,  and  corrections  involved  in  modeling 
the  ocean  floor  as  a  flat  Lagrangian  solid  bottom  in 
Dyna_N.  This  method  of  bottom  modeling  has  a  distinct 
advantage  over  Gemini  in  that  the  application  of  a 
contoured  Lagrangian  solid  bottom  model  is  has  no 
difference  from  that  of  a  flat  solid  bottom.  In  most  deep 
water  UNDEX  events  the  contour  of  the  ocean  bottom  is 
trivial  as  the  bottom  reflection  of  the  shockwave  is  of 
minimal  magnitude.  This  assumption  is  not  true  for  littoral 
waters.  The  ability  to  model  contoured  shallow  water 
environments  could  prove  vital  in  determining  the  true 
nature  of  UNDEX  effects  on  ships  operating  in  these  waters. 
With  this  in  mind,  five  different  bottom  contours  and  one 
flat  bottom  model  were  developed  and  simulated  with  the  FSP 
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serving  as  the  ship  model.  The  goal  was  to  determine  the 
effect  which  bottom  geometry  had  upon  the  response  of  both 
the  fluid  domain  and  the  FSP. 


A.  DESCRIPTION  OF  BOTTOM  SHAPES 

The  profiles  of  each  of  the  bottom  contours  used  are 
displayed  in  Figure  38.  Table  17  contains  the  values  of  the 


variables 

referenced 

rn  Figure 

38  . 

In 

all  cases  the 

standard 

charge  of  60 

pounds  of 

HBX-1 

at 

a  depth  of  ten 

meters  was  placed  with  six  meters  of  lateral  separation 
from  the  FSP.  As  discussed  in  Section  III.E,  the  dimensions 
of  the  fluid  domain  were  84  meters  in  both  the  X  and  Y 
directions  with  the  charge  placed  in  the  center.  The  depth 
of  the  fluid  extended  to  one  meter  beyond  the  lowest  point 
of  the  water-soil  interface  for  each  contoured  model.  Every 
simulation  was  run  out  to  a  full  second  to  ensure  that  the 
response  beyond  the  first  bubble  pulse  was  captured. 


Table  17.  Values  of  Contoured  Bottom  Shapes 


Figure 

XI  (m) 

X2  (m) 

Y  (m) 

Z1  (m) 

Z2  (m) 

Z3  (m) 

A 

Deep  V 

42.01 

42 .01 

84 . 02 

25.00 

20.00 

5.00 

B 

Inv .  V 

42.01 

42 .01 

84.02 

25.00 

20.00 

5.00 

C 

U  Channel 

21.01 

42 .01 

84 . 02 

10.00 

25.00 

5.00 

D 

Ramped 

84.02 

N/A 

84 . 02 

25.00 

20.00 

5.00 

E 

Anechoic 

84.00 

N/A 

84 . 00 

30.00 

10.00 

5.00 

F 

Flat 

84.02 

N/A 

84 . 02 

35.00 

5.00 

N/A 
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Figure  38.  Profile  Views  of  Contoured  Bottoms 


B.  BOTTOM  CONTOUR  EFFECT  ON  FSP  RESPONSE 

The  first  point  of  comparison  between  the  six  bottom 
types  is  the  vertical  velocity  response  of  the  closest 
point  of  the  FSP  for  the  first  800  milliseconds  in  Figure 
39.  The  velocity  response  has  been  filtered  to  eliminate 
the  all  of  the  response  data  which  had  a  frequency  greater 
than  250  Hz.  Initially  the  velocities  rise  sharply  as  the 
initial  shockwave  impacts  the  FSP.  After  approximately  50 
milliseconds  the  initial  transient  response  subsides, 
resulting  in  a  velocity  which  steadily  decreases.  Near  550 

milliseconds,  the  first  bubble  pulse  impacts  the  structure 
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causing  a  sharp  increase  in  vertical  velocity.  This  trend 
is  expected  in  nearly  all  surface  ship  responses  to  UNDEX 
events.  On  the  whole,  every  bottom  contour  causes  a  nearly 
similar  FSP  response.  The  lone  exception  is  the  anechoic 
pyramid  bottom  that  demonstrates  a  vertical  velocity 
noticeably  lower  than  the  rest. 
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Figure  39.  Velocity  Response  of  FSP  Closest  Point  to 

Charge  for  Each  Contoured  Bottom  (0-800msec) 

The  primary  concern  in  all  of  the  simulations  is  the 
effect  that  the  bottom  reflection  has  upon  the  FSP.  In  all 
six  cases,  the  bottom  reflection  should  reach  the  FSP 
within  the  first  50  milliseconds.  With  this  in  mind,  a 
closer  examination  was  conducted  of  the  early  vertical 
velocity  time-history  in  Figure  40.  The  initial  velocity 
responses  are  identical  in  all  cases  for  the  first  50 
milliseconds.  Thus  the  bottom  reflections  in  all  cases 
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appear  to  have  very  little  effect  on  the  initial  velocity 
response  of  the  FSP.  Even  after  50  milliseconds  the 
responses  show  very  little  separation. 


Series  7  -  Z  "Velocity  of  FSP 

Location:  Clossssit  Point 


T  i  me  (  ms  ) 
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Figure  40.  Velocity  Response  of  FSP  Closest  Point  to 
Charge  for  Each  Contoured  Bottom  (0-100msec) 

The  lack  of  velocity  differences  indicates  that  there 
must  be  only  minor  differences  in  the  incident  pressure  the 
FSP  feels  in  all  cases.  Before  examining  the  pressure 
history  of  the  fluid  below  the  FSP,  a  brief  understanding 
of  the  entire  pressure  history  response  is  required.  Figure 
41  shows  the  pressure  time  history  for  the  first  second  of 
the  UNDEX  simulation  and  is  representative  of  nearly  every 
UNDEX  event.  Initially,  there  is  a  large  pressure  rise  from 
the  initial  shockwave.  Assuming  a  bulk  cavitation  zone 
forms,  when  the  cavitation  zone  collapses  it  emits  a 
smaller  pressure  wave  known  as  the  cavitation  closure.  This 
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normally  occurs  within  the  first  hundred  milliseconds.  All 
of  the  subsequent  pressure  spikes  are  a  result  of  bubble 
pulses . 
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Pressure  History  Near  Surface 


Figure  41.  Example  UNDEX  Pressure  History 


The  pressure  time  histories  for  a  point  directly  below 
the  center  of  the  FSP  for  all  cases  was  plotted  in  Figure 
42  for  the  first  hundred  milliseconds  of  the  simulation. 
Indeed,  there  are  no  outliers  in  the  pressure  data.  All  of 
the  cases  followed  the  same  pressure  time  history  until 
approximately  45  milliseconds,  before  they  began  to 
diverge.  Even  after  45  milliseconds,  the  differences  were 
only  minimal  and  contained  no  spikes  in  the  pressure. 
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Sei'ies  7:  Pressure 


Location:  Center  of  Bottom 


7-01  Deep  V  Channel  7-02  Inverted  V  Climuiel 

7-03:  Deep  I.'  Channel  - 7-04  Slanted  Bottom 

7-0?:  Aneclioic  PtTiiinids  - 7-06:  Flat  Bottom 


Figure  42.  Bottom  Contour  Effect  on  Pressure  Directly 

Beneath  FSP 

While  it  appears  from  the  analysis  of  the  FSP  response 
that  the  bottom  contour  has  little  effect,  further 
investigation  of  the  fluid  domain  response  highlighted  the 
differences  and  provided  an  explanation  for  its  lack  of 
effects  on  the  FSP. 

C .  BOTTOM  CONTOUR  EFFECT  ON  FLUID  DOMAIN  AND  BUBBLE 

DYNAMICS 

The  lack  of  pressure  and  velocity  differences  directly 
below  the  FSP  between  the  different  bottom  contours  does 
not  imply  that  there  were  no  differences  in  the  fluid 
domain  response.  A  study  of  the  pressure  history  at  a  point 
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five  meters  beneath  the  FSP  was  conducted.  Figure  43  is 
focused  on  the  time  period  of  the  initial  shockwave 
propagation.  Figure  44  is  centered  on  the  first  bubble 
pulse  response. 


Figure  43.  Bottom  Contour  Effect  on  Pressure  5m  Below 

FSP  (20-50msec) 


The  two  most  similar  responses  are  the  flat  bottom  (7- 

06)  and  the  Deep  U  Channel  (7-03)  .  These  two  contours 

provided  nearly  identical  bottom  reflections  of  2.7  x  105  Pa 

at  37.5  milliseconds.  This  is  expected  since  the  point  of 

reflection  on  the  bottom  in  both  cases  was  at  the  same 

depth  and  same  zero  inclination.  The  Slanted  bottom  (7-04) 

reflection  returned  slightly  earlier  at  36  milliseconds 

with  a  greater  magnitude  of  3  x  105  Pa.  As  the  point  of 
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reflection  in  this  case  was  slightly  shallower  than  the 
flat  bottom,  the  reflection  should  be  sooner  and  of  greater 
magnitude.  This  same  principle  applies  for  the  Deep  V 
Channel  (7-01),  where  the  point  of  reflection  is  at  a 
greater  depth.  As  expected  the  reflection  was  delayed.  Its 
pressure  should  have  been  smaller  than  the  flat  bottom,  but 
was  of  equal  magnitude  instead.  The  Inverted  V  Channel  (7- 
02)  had  the  shallowest  reflection  point  and  returned  the 
reflection  the  quickest  at  26  milliseconds.  Unexpectedly, 
even  with  the  shortest  distance  to  travel,  the  Inverted  V 
Channel  gave  the  second  lowest  pressure  magnitude.  Lastly, 
the  Anechoic  Pyramid  bottom  returns  the  smallest  reflection 


wave 

,  whose  magnitude 

is 

less 

than 

hydrostatic 

pressure . 

The 

lack  of  a  bottom 

reflection  of 

significance 

explains 

why 

this  simulation 

had 

the 

lowest  vertical 

velocity 

through  the  simulation  in  Figure  39. 

The  pressure  peaks  from  45  to  50  milliseconds  were 
verified  to  be  the  result  of  cavitation  closure  with  the 
exception  of  the  Deep  V  Channel  bottom  reflection  at  44 
milliseconds.  The  pressure  peak  in  the  Deep  U  Channel 
response  at  24  milliseconds  is  the  shockwave  reflection 
from  the  vertical  face  of  the  channel.  Its  motion  was 
primarily  horizontal  and  did  not  propagate  towards  the  FSP. 

The  fluid  domain  analysis  additionally  examined  the 
differences  in  the  first  bubble  pulse  due  to  the  bottom 
contour.  The  resulting  pressure  histories  are  shown  in 
Figure  44.  While  the  magnitude  appears  unaffected,  the 
timing  of  the  pulse  shows  significant  differences  for  the 
various  contours.  Once  again,  the  two  contours  that  showed 
the  most  significant  difference  from  the  flat  bottom  were 
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the  Anechoic  Pyramids  (7-05)  and  the  Inverted  V  Channel  (7- 
02) .  These  differences  in  pulse  timing  can  have  significant 
impact  on  the  total  effect  of  a  charge  on  a  target.  The 
majority  of  UNDEX  weapons  are  designed  such  that  their 
bubble  pulses  are  timed  to  excite  the  natural  bending 
frequency  of  the  ship  thereby  causing  resonance  and 
increased  damage  [21]  .  While  the  variation  in  frequency  of 
the  bubble  response  of  a  60  pound  charge  is  small  due  to 
the  bottom  contour,  most  UNDEX  threats  are  one  or  two 
orders  of  magnitude  larger  which  have  the  potential  for  a 
wider  variation  in  pulse  frequency. 
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Figure  44.  Bottom  Contour  Effect  on  Pressure  5m  Below 

FSP  ( 550-600msec) 
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Table  18.  Contoured  Bottom  Effect  Characterization  at  5m  below  FSP 


Initial  Bottom  Reflection 

Lst  Bubble  Pulse 

Contour 
(Run  #) 

Pressure 
(%  of  Max) 

Path  Distance 
(cm) 

Pressure 

(Pa) 

Arrival  Time 
(msec) 

Frequency 

(Hz) 

Flat  Bottom 
(7-06) 

12 . 6 

7037 

1.07E+06 

566 

1 .77 

Deep  V  Channel 
(7-01) 

12 . 6 

8620 

1 . 08E  +  06 

567 

1.76 

Inverted  V  Channel 
(7-02) 

9.2 

5460 

1 . 12E  +  06 

575 

1 .74 

Deep  U  Channel 
(7-03) 

13.6 

7037 

1 . 11E  +  06 

573 

1.75 

Slanted  Bottom 
(7-04) 

14 . 6 

6835 

1  .  10E  +  06 

568 

1.76 

Anechoic  Pyramids 
(7-05) 

4 . 9 

7037 

1 . 09E  +  06 

581 

1 . 72 
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The  findings  in  this  section  have  been  summarized  in 
Table  18.  For  the  initial  bottom  reflection  the  pressure  is 
given  as  a  percentage  of  the  maximum  shockwave  pressure 
that  the  test  point  received  from  the  initial  blast.  The 
path  distance  is  the  length  of  the  bottom  reflection 
propagation  path.  The  second  section  of  the  table  includes 
the  pertinent  data  regarding  the  first  bubble  pulse. 

D.  INDIVIDUAL  BOTTOM  CONTOUR  ANALYSIS 

This  section  investigates  each  bottom  contour  response 
individually  to  determine  why  the  initial  FSP  response  is 
unaffected  by  the  bottom  contours  that  have  been  shown  to 
significantly  affect  the  fluid  domain. 

1 .  Flat  Bottom 

Examining  Figure  45,  a  clear  picture  of  the  bottom 
reflection  along  with  the  bulk  cavitation  zone  at  28 
milliseconds  is  formed.  A  black  dashed  contour  line  was 
added  to  the  figure  to  highlight  the  wave  front.  The  bottom 
reflection  should  have  impacted  the  FSP  at  approximately  40 
milliseconds.  Figure  46  shows  that  this  is  not  the  case  and 
provides  a  better  understanding  of  why  the  bottom 
reflection  has  little  effect  on  the  FSP  initially.  As  the 
bottom  reflection  travelled  vertically  through  the  water 
column,  it  impacted  the  existing  bulk  cavitation  zone. 
Recall  that  bulk  cavitation  is  created  when  a  compressive 
wave  is  incident  on  a  free  surface.  Here  the  low  pressure 
bulk  cavitation  zone  serves  as  the  free  surface  for  the 
bottom  reflection.  The  result  is  the  formation  of  a  second 
cavitation  zone  beneath  the  first.  This  indicates  that  the 
bulk  cavitation  zone  acts  as  a  buffer  for  the  FSP. 
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RunID  7-06:  Flat  Bottom 


Euler:  pressure  (Pa) 
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Figure  45. 
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RunID  7-06:  Flat  Bottom 
Euler:  pressure  (Pa) 

Lagrange:  Vertical  Stress  (dyne /cm* 2) 
Time  :  40.03000E-03 


*  C*4  291.  1.174 


Figure  46.  Pressure  Distribution  for  Run  7-06:  Flat 

Bottom  at  40  msec 
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2. 


Deep  V  Channel 


The  fluid  response  at  32  milliseconds  in  the  Deep  V 
Channel  showed  significant  differences  from  the  flat 
bottom.  In  this  case  the  contour  created  two  bottom 
reflections  seen  in  Figure  47.  At  the  convergence  of  these 
two  reflections,  a  vertically-moving,  high  pressure  zone 
formed.  Due  to  the  symmetry  of  the  bottom  contour  about  the 
charge  location,  this  high  pressure  zone  collided  with  the 
gas  bubble  in  Figure  48,  which  dissipated  its  energy,  and 
did  not  immediately  affect  the  FSP.  In  this  case  there  is 
no  evidence  of  a  second  bulk  cavitation  zone  forming. 


RunID  7-01:  Deep  V  Channel 
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Figure  47.  Pressure  Distribution  for  Run  7-01:  Deep  V 

Channel  at  32  msec 
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RunlD  7-01:  Deep  V  Channel 


Euler:  Pressure  (Pa) 

Lagrange:  Vertical  Stress  (dyne/cmA2) 
Time  :  41.05000E-03 


-5  00006*0$ 
-i  oooot-oe 
-1  5GOOC*C6 
-7  00006*06 

-?$oooe*o» 

-300006*06 
-4  0000E-06 
■500006*06 
•1  00006*07 

•i  soooe*o7 

•7  0000E*C7 

Mn  -8  97096*06 
tfNooe  19101 
U*>  100*36*06 
«N0<M  3914® 


Pt**) 


B 


1  S000E»07 
1  00006*07 
500006*06 
.'50006*06 
6  00006*05 
$00006*05 
4  00006*05 
300006*0$ 
700006*0$ 
1  01006*06 
990006*04 
$70006*03 
000006*00 


Mn  $00006*03 
tf  C«4  177.  1. 
Ms*  6  $3176*0$ 
MC4®  10.  1.173 


Figure  4  8 . 


Pressure  Distribution  for  Run  7-01: 
Channel  at  41  msec 


Deep  V 


3. 


Inverted  V  Channel 


Unlike  the  Deep  V  Channel  where  high  pressure 
convergence  zones  were  created,  the  Inverted  V  Channel 
causes  the  bottom  reflection  in  Figure  49  to  spread  after 
only  20  milliseconds.  The  dispersion  of  the  bottom 
reflection  helps  explain  the  reduced  pressure  noted  in 
Figure  43.  Once  again,  a  new  bulk  cavitation  zone  is 
observed  after  the  bottom  reflection  is  incident  on  the 
original  cavitation  zone  in  Figure  50. 
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Bulk  Cavitation 


RunID  7-02:  Inverted  V  Channel 
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Figure  49. 
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RunID  7-02:  Inverted  V  Channel 


Euler:  Pressure  (Pa) 

Lagrange:  Vertical  Stress  (dyne/cmA2) 
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Pressure  Distribution  for  Run  7-02:  Inverted 
V  Channel  at  28  msec 
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Figure  50 


4. 


Deep  U  Channel 


The  pressure  contour  plot  for  the  Deep  U  Channel  at  25 
milliseconds  is  provided  in  Figure  51.  The  view  of  the 
bottom  reflection  is  ambiguous  due  to  the  presence  of 
reflections  from  the  vertical  sidewalls.  Unlike  the 
previous  three  cases  where  the  bulk  cavitation  zone  acted 
as  a  buffer,  it  appears  in  Figure  52,  that  the  cavitation 
zone  disappeared  just  prior  to  or  as  a  result  of  the  bottom 
reflection.  This  allowed  an  increase  in  pressure,  shown  in 
green,  to  be  seen  just  below  the  hull  of  the  FSP.  Though  it 
is  present,  this  pressure  is  less  than  twice  the  value  of 
hydrostatic  pressure. 
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Figure  51.  Pressure  Distribution  for  Run  7-03:  Deep  U 

Channel  at  25msec 
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Figure  52.  Pressure  Distribution  for  Run  7-03:  Deep  U 

Channel  at  39  msec 


5 .  Slanted  Bottom 

Figure  53  is  the  pressure  contour  response  for  the 
Slanted  Bottom  simulation  at  28  milliseconds.  The  bottom 
reflection  looks  similar  to  that  of  the  flat  bottom  with  no 
areas  of  convergence  or  dispersion.  In  Figure  54,  the 
formation  of  a  second  bulk  cavitation  zone  indicates  that 
the  bottom  reflection  does  not  impact  the  FSP. 
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RunID  7-04:  Slanted  Bottom 
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Figure  53. 


Pressure  Distribution  for  Run  7-04:  Slanted 
Bottom  at  28  msec 


RunID  7-04:  Slanted  Bottom 
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Figure  54.  Pressure  Distribution  for  Run  7-04:  Slanted 

Bottom  at  37  msec 
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6. 


Anechoic  Pyramids 


The  purpose  of  an  anechoic  surface  is  to  minimize  or 
eliminate  reflections.  The  use  of  a  bottom  contour  of 
anechoic  pyramids  accomplished  this  function.  Figure  55 
highlights  that  the  bottom  reflections  created  by  the 
anechoic  surface  were  weak  and  scattered.  By  40 
milliseconds  in  Figure  56,  there  are  no  distinct  pressure 
waves  remaining  in  the  fluid  domain.  This  corresponds  well 
to  the  minimal  pressure  wave  response  that  was  noted  in 
Figure  43. 


RunID  7-05:  Anechoic  Pyramids 
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RunID  7-05:  Anechoic  Pyramids 


Euler:  Pressure  (Pa) 

Lagrange:  Vertical  Stress  (dyne/cmA2) 
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Figure  56.  Pressure  Distribution  for  Run  7-05:  Anechoic 

Pyramids  at  40  msec 


The  anechoic  surface  demonstrates  a  unique  fluid 
response  where  pockets  of  cavitation  form  at  the  water-soil 
interface.  A  review  of  the  animated  pressure  contour 
response  shows  that  over  time  these  cavitation  pockets 
migrate  from  the  bottom  upward  around  the  bubble.  These  low 
pressure  cavitation  zones  do  not  slow  the  bubble  expansion 
at  the  same  rate  as  the  normal  fluid  domain.  This  allows 
the  bubble  to  expand  over  a  longer  time,  thereby  delaying 
the  first  pulse. 

E.  BULK  CAVITATION  AS  A  BUFFER  ZONE 

It  is  now  clear  that  the  lack  of  initial  FSP  response 
to  the  different  bottom  contours  was  due  to  the  buffer 
provided  by  the  bulk  cavitation  zone.  This  effect  is 
strictly  a  by-product  of  this  particular  simulation 
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geometry,  upon  which  every  simulation  was  based.  The  charge 
weight,  depth,  target,  and  lateral  separation  between  the 
two  were  chosen  to  mimic  the  MIL-S-901D  testing  that 
utilizes  the  FSP  as  a  shock  test  platform  for  shipboard 
equipment  [5]  .  The  chosen  bottom  depth  was  an  average  of 
several  shallow  water  environments. 

Several  factors  can  be  varied  in  order  to  diminish  the 
buffer  effect  of  the  bulk  cavitation  zone.  As  noted  in 
Figure  43,  the  cavitation  closure  for  this  particular 
charge  occurs  at  approximately  50  milliseconds.  An  increase 
in  bottom  depth  would  allow  the  bottom  reflection  to  return 
to  the  surface  after  the  cavitation  zone  closed.  Thus  the 
bottom  reflection  would  reload  the  structure.  The  draft  of 
the  FSP  is  only  1.2  meters.  If  the  Littoral  Combat  Ship 
(LCS)  with  a  draft  of  4.5  meters  were  used  the  keel  would 
be  below  the  buffer  zone.  A  lateral  shift  of  the  target  to 
a  point  outside  the  bulk  cavitation  zone,  it  would  feel  the 
full  effect  of  the  bottom  reflection  no  matter  the  ocean 
depth.  As  charge  size  and  depth  increase,  the  lateral 
extent  of  the  bulk  cavitation  zone  does  as  well  to  a  point. 
Lastly,  while  the  cavitation  zone  acted  as  a  buffer  for  the 
bottom  reflection  of  a  60  pound  charge,  a  larger  charge 
with  a  larger  bottom  reflection  might  be  able  to  collapse 
the  bulk  cavitation  zone  and  still  be  energetic  enough  to 
significantly  impact  the  FSP. 
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VII.  FINAL  REMARKS 


A.  CONCLUSIONS 

A  method  to  model  the  ocean  bottom  as  a  Lagrangian 
solid  was  developed  for  comparison  to  the  current  Euler 
fluid  bottom  modeling  approaches.  Several  sets  of  boundary 
and  initial  conditions  were  simulated  to  determine  which 
combination  introduced  the  least  artificiality  to  the  fluid 
domain  solution.  If  an  elastic  or  elasto-plastic  model  was 
used,  it  was  then  possible  to  apply  non-reflecting  boundary 
segments  and  nodal  reaction  forces  at  the  boundaries  along 
with  Dyna_N  Prestress.  The  non-reflecting  boundary  segments 
allowed  the  Lagrangian  solid  bottom  to  act  as  a  semi¬ 
infinite  domain,  thereby  eliminating  retransmission  waves. 
The  Dyna_N  Prestress  imposed  the  hydrostatic  loading  and 
deformation  on  the  bottom  in  order  to  minimize  the 
magnitude  of  the  initial  bottom  wave.  This  combination 
provided  an  accurate  and  efficient  solution.  However,  the 
non-reflecting  boundary  segments  and  Dyna_N  Prestress  have 
limitations  that  prohibit  their  application  to  more 
complicated  soil  material  models. 

The  validation  of  the  Lagrangian  solid  bottom  model 
was  completed  by  comparing  its  fluid  domain  response  to  two 
existing  bottom  modeling  methods.  The  first  bottom  modeling 
method  was  a  purely  reflective  Eulerian  boundary  condition. 
The  second  method  was  the  use  of  an  Eulerian  equation  of 
state  for  a  generic  clay  soil.  The  simulation  geometry  was 
consistent  with  the  MIL-S-901D  shock  testing  utilizing  the 
Floating  Shock  Platform.  Although  the  magnitude  was 
greater,  the  bottom  reflection  of  the  Lagrangian  solid 
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bottom  model  had  the  same  characteristic  shape  of  the  clay 
soil  reflection.  The  first  bubble  pulse  occurred  13 
milliseconds  sooner  with  the  solid  bottom  than  with  the 
Euler  clay  soil.  Overall,  the  fluid  response  of  the 
simulation  with  a  Lagrangian  solid  bottom  was  consistent 
with  the  Euler  bottom  modeling  method. 

One  of  the  distinct  advantages  of  using  a  Lagrangian 
solid  bottom  is  the  ability  to  model  contoured  bottom 
shapes.  Six  contoured  solid  bottom  models  were  developed  to 
investigate  bottom  contour  effects  on  shallow  water  UNDEX 
events.  The  initial  analysis  of  the  FSP  response  showed 
only  slight  differences  between  the  various  contour  models. 
This  was  caused  by  the  buffer  created  by  the  bulk 
cavitation  zone.  The  effect  was  specific  only  to  the 
particular  geometry  selection.  Modifications  of  the  charge 
size,  target  separation,  or  bottom  depth  could  diminish  the 
effect,  but  it  was  left  for  future  study.  Further 
investigation  of  the  fluid  domain  response  revealed  that 
there  were  indeed  significant  differences  between  the 
initial  bottom  reflections  for  the  different  contours.  The 
most  important  bottom  contour  effect  was  the  distortion  to 
the  gas  bubble  and  its  associated  first  pulse  timing.  These 
changes  could  have  severe  implications  in  the  case  of 
undersea  weapons  designed  to  take  advantage  of  ship 
whipping . 

B .  FURTHER  RESEARCH 

While  the  preprocess  required  to  implement  the 
Lagrangian  solid  bottom  model  can  be  time-intensive,  the 
potential  benefits  are  considerable. 
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Application  of  Dyna_N  material  model  types  16,  45,  and 
65,  which  were  created  to  model  geologic  material,  to  the 
Lagrangian  solid  bottom  should  be  investigated.  A  more 
accurate  soil  model  could  further  improve  the  accuracy  of 
the  fluid  domain  response. 

A  recent  DYSMAS  code  revision,  which  was  unavailable 
during  this  research,  allows  for  the  use  of  viscous  fluids. 
This  viscous  code  was  developed  in  order  to  more  accurately 
model  the  ocean  bottoms.  A  comparison  of  capabilities  could 
be  conducted  between  this  new  code  and  the  Lagrangian  solid 
bottom  model,  especially  in  the  case  of  charges  placed 
close  to  the  sea  floor  where  the  effects  of  cratering  are 
highly  likely. 

Although  it  was  designed  for  application  in  shallow 
water  simulations,  the  Lagrangian  solid  bottom  model  has 
potential  applications  in  deep  water  simulations. 
Currently,  all  fluid  domain  simulations  in  which  gravity  is 
present  require  the  application  of  a  reflective  bottom 
boundary  condition.  In  deep  water  simulations,  a  sufficient 
fluid  depth  is  added  to  the  domain  in  order  to  minimize  the 
effect  of  this  reflective  boundary.  Recalling  that  the 
behavior  of  waves  at  an  interface  is  determined  by  the 
ratio  of  the  product  of  each  materials  density  and  sound 
speed,  it  could  be  possible  with  an  elastic  material  model 
to  create  an  interface  ratio  of  one.  This  would  allow  the 
wave  energy  to  be  completely  absorbed  into  the  Lagrangian 
solid  bottom.  The  use  of  NRB  segments  would  then  allow  for 
the  dissipation  of  the  wave.  Thus  the  depth  of  the  fluid 
could  be  decreased  to  only  the  area  of  interest  around  the 
gas  bubble  and  target.  If  the  material  properties  could  not 
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be  matched,  then  creating  a  contoured. 


anechoic  bottom 


surface 

provide 


could  also  eliminate  the  bottom  reflection  and 
a  more  accurate  deep  water  solution. 
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APPENDIX.  INDEX  OF  SIMULATIONS 


The  following  table  is  an  index  of  the  UNDEX 
simulations  that  were  conducted  in  the  course  of  this 
research.  The  pertinent  data  for  the  charge,  fluid 
geometry,  FSP  position,  solid  bottom  dimensions,  and  solid 
bottom  boundary  and  initial  conditions  are  listed  for  each 
simulation.  For  additional  simulation  input  data  and 
results,  contact  the  Shock  and  Vibration  Computational 
Laboratory  at  the  Naval  Postgraduate  School. 
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Run  ID# 

1-00 

1-01 

3-01 

3-02 

3-03 

3-07 

3-08 

4-01 

4-02 

4-04 

4-07 

5-03 

Description 

Wall  Boundary  Condition 

for  Ocean  Bottom 

Gemini  Soil  EOS  for 

Ocean  Bottom 

Long-Time  Steady-State 

Simulation  (2.5  secs) 

Long-Time  Steady-State 

Simulation  w/  Ramped 

Load 

DR  of  Structural  Bottom 

DR  of  Structural  Bottom 

and  FSP 

Dyna_N  Prestress  of 

Structural  Bottom 

Retransmission  of 

Nodally  Constrained 

Bottom  5m 

Retransmission  of 

Nodally  Constrained 

Bottom  10m 

Effect  of  N  RB  Segments 

on  Retransmission 

Effect  of  NRBs  and  Nodal 

Forces  on 

Retransmission 

Bottom  Contour  in  Euler 

(2D)  Proof  of  Concept 

Simulation  Time  (sec) 

1.000 

1.000 

2.500 

2.500 

N/A 

N/A 

0.100 

0.100 

0.100 

0.100 

0.100 

0.100 

Charge 

Mass  (kg) 

27.2 

27.2 

N/A 

N/A 

N/A 

N/A 

27.2 

27.2 

27.2 

27.2 

27.2 

N/A 

Depth  (m) 

10.0 

10.0 

N/A 

N/A 

N/A 

N/A 

10.0 

10.0 

10.0 

10.0 

10.0 

N/A 

X  (m) 

84.0 

84.0 

84.0 

84.0 

N/A 

N/A 

84.0 

84.0 

84.0 

84.0 

84.0 

10.0 

Y(m) 

84.0 

84.0 

84.0 

84.0 

N/A 

N/A 

84.0 

84.0 

84.0 

84.0 

84.0 

N/A 

Fluid  Domain  Extent 

Z(m) 

45.0 

50.0 

46.0 

46.0 

N/A 

N/A 

46.0 

46.0 

46.0 

46.0 

46.0 

10.0 

Max.  Water  Depth  + 

35.0  + 

35.0  + 

35.0  + 

35.0  + 

35.0  + 

35.0  + 

35.0  + 

35.0  + 

35.0  + 

35.0  + 

35.0  + 

8.0  + 

Overlap  (m) 

0.0 

5.0* ** 

1.0 

1.0 

0.0 

0.0 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0* 

X  (m) 

853.4 

853.4 

N/A 

N/A 

N/A 

853.4 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

FSP  Position 

Y(m) 

0.0 

0.0 

N/A 

N/A 

N/A 

0.0 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

Z(m) 

-121.9 

-121.9 

N/A 

N/A 

N/A 

-121.9 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

X  (m) 

N/A 

N/A 

84.00 

84.00 

84.00 

84.00 

84.00 

84.00 

84.00 

84.00 

84.00 

N/A 

Solid  Bottom  Size 

Y(m) 

N/A 

N/A 

84.00 

84.00 

84.00 

84.00 

84.00 

84.00 

84.00 

84.00 

84.00 

N/A 

Z(m) 

N/A 

N/A 

5.00 

5.00 

5.00 

5.00 

5.00 

5.00 

10.00 

5.00 

5.00 

N/A 

Interface  Segments 

N/A 

N/A 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

N/A 

Solid  Bottom 

Boundary 

Conditions 

Fixed  Nodes 

N/A 

N/A 

Yes, 

X/Y/Zmin 

Yes, 

X/Y/Zmin 

No 

Yes, 

X/Y/Zmin 

No 

Yes, 

X/Y/Zmin 

Yes, 

X/Y/Zmin 

No 

No 

N/A 

Reaction  Forces 

N/A 

N/A 

No 

No 

Yes, 

ANSYS 

No 

Yes, 

Manual 

No 

No 

No 

Yes, 

Manual 

N/A 

NRB  Segments 

N/A 

N/A 

Yes 

Yes 

No 

No 

Yes 

No 

No 

Yes 

Yes 

N/A 

Counter  Back  Pressure 

N/A 

N/A 

Yes 

Yes 

Yes 

No 

No 

No 

No 

No 

No 

N/A 

Solid  Bottom  Initial 

Conditions 

Dyna_N  Prestress 

N/A 

N/A 

No 

No 

No 

No 

Yes, 

Manual 

No 

No 

No 

No 

N/A 

*  indicate  overlap  was  actual  euler  clay  material 

**  indicates  average  depth,  See  Table  17  and  Figure  38  for  more  accurate  description 
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Run  ID# 

7-01 

7-02 

7-03 

7-04 

7-05 

7-06 

8-01 

8-02 

8-03 

8-04 

8-05 

8-06 

Description 

Contoured  Structural 

Bottom  -  Deep  V 

Contoured  Structural 

Bottom  -  Inverted  V 

Contoured  Structural 

Bottom  -  Deep  U 

Channel 

Contoured  Structural 

Bottom  -  Slanted 

Contoured  Structural 

Bottom  -  Anechoic 

Pyramids 

Flat  Bottom 

Effect  of  Bottom 

Thickness  w/  Prestress 

2m 

Effect  of  Bottom 

Thickness  w/  Prestress 

5m 

Effect  of  Bottom 

Thickness  w/  Prestress 

10m 

Effect  of  Bottom 

Thickness  w/o  Prestress 

2m 

Effect  of  Bottom 

Thickness  w/o  Prestress 

5m 

Effect  of  Bottom 

Thickness  w/o  Prestress 

10m 

Simulation  Time  (sec) 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

0.100 

0.100 

0.100 

0.100 

0.100 

0.100 

Charge 

Mass  (kg) 

27.2 

27.2 

27.2 

27.2 

27.2 

27.2 

27.2 

27.2 

27.2 

27.2 

27.2 

27.2 

Depth  (m) 

10.0 

10.0 

10.0 

10.0 

10.0 

10.0 

10.0 

10.0 

10.0 

10.0 

10.0 

10.0 

Fluid  Domain  Extent 

X  (m) 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

Y  (m) 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

84.0 

Z(m) 

56.0 

56.0 

46.0 

46.0 

41.0 

46.0 

46.0 

46.0 

46.0 

46.0 

46.0 

46.0 

Max.  Water  Depth  + 
Overlap  (m) 

45.0  + 

1.0 

45.0  + 

1.0 

35.0  + 

1.0 

35.0  + 

1.0 

40.0  + 

1.0 

35.0  + 

1.0 

35.0  + 

1.0 

35.0  + 

1.0 

35.0  + 

1.0 

35.0  + 

1.0 

35.0  + 

1.0 

35.0  + 

1.0 

FSP  Position 

X  (m) 

853.4 

853.4 

853.4 

853.4 

853.4 

853.4 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

Y(m) 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

Z(m) 

-121.9 

-121.9 

-121.9 

-121.9 

-121.9 

-121.9 

N/A 

N/A 

N/A 

N/A 

N/A 

N/A 

Solid  Bottom  Size 

X  (m) 

84.02 

84.02 

84.02 

84.02 

84.00 

84.02 

84.02 

84.02 

84.02 

84.02 

84.02 

84.02 

Y(m) 

84.02 

84.02 

84.02 

84.02 

84.00 

84.02 

84.02 

84.02 

84.02 

84.02 

84.02 

84.02 

Z(m) 

**  15.00 

**  15.00 

**  15.00 

**  10.00 

10.00 

5.00 

2.00 

5.00 

10.00 

2.00 

5.00 

10.00 

Solid  Bottom 
Boundary 
Conditions 

Interface  Segments 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

Yes, 

X/Y/Zmax 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

Yes,  Zmax 

Fixed  Nodes 

No 

No 

No 

No 

No 

No 

No 

No 

No 

No 

No 

No 

Reaction  Forces 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

NRB  Segments 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Counter  Back  Pressure 

No 

No 

No 

No 

No 

No 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

Solid  Bottom  Initial 

Conditions 

Dyna_N  Prestress 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

Yes, 

ANSYS 

No 

No 

No 

*  indicate  overlap  was  actual  euler  clay  material 

**  indicates  average  depth,  See  Table  17  and  Figure  38for  more  accurate  description 
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